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EDITOR'S NOTE 


The papers presented herein have been derived primarily from speakers' summaries of talks 
presented at the Flight Mechanics/Estimation Theory Symposium held October 18 and 19, 
1977 at Goddard Space Right Center. For the sake of completeness, abstracts are included 
of those talks for which summaries were unavailable at press time. Papers included in this 
document are presented basically as received from the authors with a minimum of editing. 
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COMPRESSION OF EPHEMERIOES BY DISCRETE CHEBYSHEV APPROXIMATIONS 

H.H. Pickard 

U.S. Naval Research Laboratory 

A. Oeprit and W. Poplarchek 
University of Cincinnati 


There has been a lot of Interest lately In representing the ephenerldes of 
satellites and planets In terms of truncated polynomial series. This paper dis- 
cusses the use of Chebyshev series for this purpose and specifically a Fortran 
package which has been developed for fitting satellite orbits. The .features 
desired In any approximation are 1) the ability to compress a satellite ephemerls, 
2) the ability to represent a satellite ephemerls over several orbits, 3) 
guaranteed accuracy to within prescribed tolerance over the time Interval of 
consideration, and 4) fast processing. These features are Imposed with an eye 
towards adapting the approximation for use on microprecessor applications In 
which storage Is limited and real time processing Is required. GPS, for 
example, will require that the representation be usable not only on spacecraft 
but also by users on ships, aircraft, or portable land units. The use of a 
polynomial approximation ensures the fast processing of requirement (4). above 
since only multiplications additions, and subtractions are Involved In the 
processing. The way In which polynomial approximations can be made to satisfy 
the other three requirements above Is the subject of this paper, 

Corlo [1] demonstrates the ability of a Chebyshev polynomial series to 
represent a satellite orbit by fitting 25 points with a 24th degree poly- 
nomial. The 25 points are equally spaced, and since the degree of the 
polynomial Is exactly one less than the number of points, the polynomial 
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ffltcrpoUtes the 25 points. The X coordinate of a geosyr^hronous satellite 
over two periods of Its orbit Is shown In Fig. (1). taken from Corlo's paper. 
The corresponding error curve, which Is the difference between the true value 
of X and the polynomial approximation to X. Is plotted In Fig. (2), also 
from Corlo's paper. It Is apparent that the «rror curve Is not at all uniform. 
The approximation fits the 25 data points exactly, as It must, but the Gibbs 
effect Is striking as we see that errors of TO km. occur at the end points of 
the Interval. The desired one n»ter accuracy can only be achieved over the 
middle 24 hours of the Interval. 

This approximation can be Improved greatly by choosing at unequal 
Intervals the reference points at which we Interpolate. We can Improve even 
further (without Increasing the degree of the polynomial) by Increasing the 
number of reference points. Naturally we must abandon polynomial 1nterpola> 
tion to do this, and must use other methods such as least squares (L.S.) 
approximation or, as will be done here, linear programming (L.P.) techniques. 
The key to the whole pror«dure Is to use a non-uniform distribution of 
reference points. This causes problems, since numerical Integrators 
generally tabulate ephemerldes at equal Intervals. Therefore, the authors 
present a method which has been developed to remove this problem. 

Best Approximations 

Chevyshev conjectured that a best polynomial approximation of degree N 

to a function y exists. 

Let P(c,t) = T. cJJt) approximate y(t) 
o<3<N ^ 

£ Is the vector of coefficients Cj. 

Tj Is the Chebyshev polynomial of degree j. 

T Is linearly related to t (a £ t £ b and -1 £ t <. 1). 




Figure 1. X-Coordinate: Geosyncronous Satellite 



Figure 2. Corresponding Error Curve: 

Degree = 24 

25 Reference Points » Uniformly Spaced 
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N Is the degree of P(c»t) 

Define e(c,t) ■ y(t) - P(c^t) • 1 * i 

Find the set denoted c* such that 
Max I e(e*.t)| < Hax | e(c^t)| 
for every c and all t In a ^ t <, b. 

P* (t) • H c? Tf (t) Is then the "best approximation" of degree 
o<4<N'’ '* 

N toy(t). 

P*(t) Is characterized by the following: 

There exist at least N-»'2 points (y(t|^)»t|^) such that 
a ^ tg ^ t^ ^ ^ t|^ ^ t|^ ^ 1 K b 

e(cl» V - (-1)*' X , 

I e(c*,t) I < X 

X >.0 

The points t|^ are called "critical points". Conversely* a polynomial exhibiting 
these properties is the "best approximation" of degree N to y(t). The authors 
of [2] have developed an algorithm to move from an Initial set of reference 
points to the critical set In a finite number of Iterations. They have pro- 
grammed the algorithm In PLl and have used It In planetary applications. 
Unfortunately, the procedure requires an analytical or semi-analytical orbital 
theory, and cannot be used with the tabulated ephemerldes coming from numerical 
Integration techniques. The reason for this Is that the method requires first 
derivatives of the function being approximated, and numerical differentiation 
cannot adequately evaluate these derivatives. 


Discrete Approximations 


The polynomial laterpmatlon scheme of Corlo (In which the 






rtference points are equally spaced) does not yield a uniform approxima- 
tion of satellite orbit coordinates, though a polynomial Interpolation 
can yield good results If the reference points are selected carefully. 

Further, a "best" polynomial of given degree to a given function does exist, 
but it Is not possible to find the "best" approximation when all that Is 
known are discrete values of the function. Therefore, try an approximation 
based on fitting a number of data points which Is considerably 
larger than the desired degree of the approximating polynomial. This suggests 
that the function be approximated using a L.S. fit to the tabulated data. 
Another approach Is to fit the data using a minmax approach with an efficient 
L.P. algorithm developed by Barrodale and Phillips at the University of 
Manitoba [3]. Both L.P. and L.S. methods will yield uniform approximations. 

If the reference points are properly chosen. The Important advantage of 
using Barrodale's L.P. algorithm. Is that It automatically gives the 
maximum error In fitting the given points, which In turn gives an estimate 
of the maximum error over the entire Interval . 

The discrete approximation problem Is stated as follows: 

Let P(c,t) » 2 c^ Tj (t) approximate y(t) as In the continuous case. 

o<j<M '* •' 

Define e(£,t) as before. 

For M discrete points (y{tj,),t|^) In the Interval a £ t £ b, find a set of 
coefficients denoted c* such that 
Max |e(c*, tj^)| < Max |e(c, tj,)| 
for all £ and k » 1,2,...M. 

P*(t) ■ ^ c.*T, (t) Is then the best approximation of degree N to the 

05 l<N ^ 

M points (y(t|j),tj^). It should be noted that Barrodale's algorithm does not 
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' riqutrc « representation fn terms of Chebysbev polynomials. In general, the 

Tj above can be replaced by any set of real valued functions Oj. Barrodale's 
algorithm is a modification of the Simplex procedure for solving linear 
programing problems. The Fortran implementation of the algorithm is given 
In the ACM Transactions on Mathematical roftware£3]. A detailed description 
of Simplex metiiod used 1$ given £4} and can be obtained by writing directly 
to Barrodale or Phillips. 

Before describing the pacicage that has been developed and results that 

have been obtained with it. Figs. (3) and (4) are presented to emphasize the 

Importance of carefully choosing the reference points. The dashed curves 

in Fig. (3) Indicate the error obtained in fitting the number of .joints given 

by the abscissa with a 24th degree polynomial. The solid curves give the 

error evaluated at 500 points over the interval and are used as the ../e 

measure of the error over the continuum of the Interval. Curves 1 correspond 

to a non-uniform distribution of reference points, and Curves 2 c'*rrespond tc 

a uniform distribution. (These results are for a U hour satellite; ti>: 

error is In meters and is logarithmically scaled.) These curves show that 

a uniform distribution of points gives vury poor results when there are only 

a few reference points, while a non-uniform distribution gives reasonably 

good results. Even for-a large number of reference points, the uniform 

case never does as well as the non-uniform case, even though the predicted 

error Is always less for the uniform case. In 'order to get the lowest degree 

which will give the desired accuracy, with a minimum number of reference 

points, and with a rellaole estimate of the true error ever the Interval, it 

Is 1ng>erat1ve that the reference points be non-uni formly spaced. 

» 

The error curves of Fig. (4) also evaluate the performance of polynomial 
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Figure 4, Error Curves for 24th Degree Fit 








approximations under different situations. The fit is for two periods 
of a satellite; the error is in meters and is logarithmically scaled. 

The first curves in each set correspond to a non-uniform distribution. 

The approximation for 25 reference points is seen to be reasonably good. 

That for 60 points is close to the best approximation* since the error 
curve ripples uniformly over the interval. The second curves in each 
set correspond to a uniform distribution. For 25 points, the error is 
excessively large over a large portion of the interval. For 60 points, 
the error is much better over much of the interval, but still rises 
sharply near the endpoints. The third figure in each set is for a non- 
uniform distribution, but L.S. fitting is used rather than the L.P. 
method. The L.S- fit is certainly not any better than the L.P. 
fit, and thus, the computational superiority of the L.P. method makes 
it the method of choice. 

The NRL Package that has been built around this algorithm is tailored 
to the problem of fitting satellite orbital elements. The main input to 
the package is an ephemeris file giving time, distance, latitude, and 
longitude. The time intervals in the file may, but need not be, equal;* 

The program will automatically interpolate (using Lagrange interpolation) 
the data to get values of the elements at the desired times. The total 
time interval of consideration, the number of reference points to use, 
the desired accuracy of the interpolation, and the desired degree or range 
of degree of the approximating polynomial are entered via a separate control 
file. As output, the user receiv^'S the coefficients to construct the 
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approximating polynomial and the maximum error In fitting the reference 
points. 

Often the user will know the accuracy he wants but will not know the 
degree needed to achieve this accuracy. ■ In this case. It is possible to 
specify a range In the degree and the program will automatically Increment 
the degree until the necessary degree Is found. When running the program 
In this mode. It Is often possible for the program to determine half way 
through the fitting procedure that the specified accuracy will not be 
attained, and abort the procedure. The program then Increments the degree 
and tries again. The program also has the ability to Increment the 
degree by more than one If It appears that the accuracy attained by the 
current degree will be. much less than that desired. This method works 
well for orbits with small eccentricity, since the program converges 
rather quickly in this case. For higher eccentricities (e^O.5), there 
Is a problem. A plateau Is reached, wherein It takes large Increments In 
the degree to Improve the accuracy of the fit. 

It Is the built-in ability to estimate the error In the approximation 
which makes Barrodale's algorithm much more convenient to use than L.S. 
fitting. When L.S. fitting Is used, one obtains only the coefficients 
needed to construct the approximating polynomial. If the desired degree 
is already known, this may be sufficient, but usually some Idea of the 
accuracy of the fit Is required. With L.S., this estimate of accuracy 
must be obtained apart from the fitting procedure. From the standpoint 
of computer use, this Is awkward and Inefficient. Barrodale's algorithm. 
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on th(> other hand, gives not only the error In the fit, but a reliable 
estl'^ite of this error half way through the fitting procedure. Then, 

If ( le reference points have been properly chosen, this error In fitting 
M points provides a reliable estimate of the error over the continuum of 
the time Interval. Thus, if only the desired accuracy is known, the 
fitting routine can automatically Increment the degree until this accuracy 
Is 6 :h1eved. 

The general flow of the NRL Package Is Illustrated in Figure (5). 

The ephemerls and control files supply Input to the Interpolation block, 
which In turn prepares the reference points to be fitted. Normally, the 
ephemerls file will be much larger than the final reference set to be 
fitted. To fit all the points of the ephemerls file would be extremely 
:ostly and unnecessary. Far fewer points can usually be used If they 
are non-uniformly spaced. The interpolation block of the package obtains 
hese points quickly and accurately. The Barrodale fitting routine is 
S iow:i in the dashed box. It has been modified to allow early exit under 
certain cond1t1_..s, and a looping structure has been built around it to 
automatically increment the degree of the approximating polynomial. 

€RRW\X is a FORTRAN variable, which indicates the maximum error desired. 
ERRfiAX and the desired range of degree are both supplied through the 
control file. RESMAX is a FORTRAN variable which on termination of the 
fitting procedure ,s the maximum error of the fit. At the point of the 
first test ‘thin the dashed box (this comes about half way through the 
fitting algorithm), RESMAX is very close to but less than this maximum 
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Figure 5, Flowchart of the NRL Ephemeris Compression Package 
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error. ThuSt at this point, RESKAX Is a good Indicator of whether this 
naxlmum error will ultimately be achieved. The NRL package has used this 
feature of Barrodale's algorithm to advantage In the looping structure. 
Most of the looping can take place quickly, so that the second test (just 
outside the dashed box) Is usually satisfied when It Is reached. If so, 
the fitting procedure Is finished. If not, the looping continues auto- 
matically. 

The dasheo box In the lower right hand corner of figure (5) Is 
Included to demonstrate the procedure If the L.P. algorithm Is replaced 
with a L.S. algorithm. Neither of the two tests In the L.P. method can 
be used In a L. S. method. The entire procedure must be finished and 
the approximation evaluated at each of the points of the reference 
file to obtain the error of the fit. In the NRL package, this box 
has been Included as an option to evaluate the performance of the 
polynomial approximation. The maximum error over the entire ephemerls 
file Is calculated for comparison with RESfWX coming out of Barrodale's 
algorithm. Experience has shown that for satellites of small to 
moderate eccentricity, RESMAX Is accurate to within 10% - 20X. Thus, 
this final check can be bypassed In most cases. 

Representative examples of the use of the package are shown in Tables 
(1), (2), and (3). Table (1) shows the results for the compression of 
the 24^ satellite SMS-B. These results show that one period or less of 
a satellite can be easily fitted with a polynomial of low degree. Equally 
Important Is the fact that several orbits can be fitted with a polynomial 



Table 1. Compression of Ephemerldes for tbe 24h Satellite SMS-B 



RANGE 

DEGREE 

(METERS) 

(SEC OF ARC) (SEC OF ARC) 


0^,5 

8 

2.28 

.003 

.002 



12 

0.83 

.004 

CM 

O 

O 

* 


2<1. 

20 

3.47 

.02 

.003 


3<1. 

28 

6.17 

.04 

.004 


4<1. 

36 

9.68 

.07 

.01 



40 

33.45 

.21 

.02 



Table 2. 

Compression of Lunar Ephemeris 



aR 

dU 

4V, 

ta AS 

RANGE DEGREE 

(METERS) 

(SEC OF ARC) 

(SEC OF ARC) 

(SEC OF ARC) (SEC OF ARC) 

7'* 

7 

3.4 

.02 

.01 

1.00 1.17 

lA'* 

10 

3.70 

.05 

.02 

.13 .12 

21 d 

18 

12,70 

.01 

.01 

.08 .03 


40 

a.3 

.01 

.01 

,01 .01 


50 

9.87 

.02 

.01 

o 

o 


Table 3, Difference in Degree Needed to Fit Elliptic and Real Orbits 

1 km 

0 ( 4) 

6 { 6) 

18 (18) 

R COMPONENT 

l^m 

0 (10) 

12 (12) 

34 (37) 

1 PERIOD 

1 km 

11 (11) 

11 (11) 

25 (25) 

X CCVPOMEHT 

1 m 

15 (15) 

15 (15) 

35 (43' 

1 PERIOD 

1 km 

0 ( 8) 

12 (12) 

* 

R COMPONENT 

1 m 

15 (15) 

22 (22) 

it 

2 PERIODS 

1 km 

18 (18) 

20 (20) 

it 

X COMPONENT 

1 m 


30 (30) 

* 

2 PERIODS 
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of moderate degree. In this case, five periods of the satellite are 
approximated by a polynomial of degree 40. Table (2) gives results for 
the compression of the lunar ephemeris. This table again demonstrates 
the ability of polynomial approximation to represent more than one 
period of an orbit. A polynomial of degree SO can represent 2 periods 
of the moon with an error of less than 10 meters, which is better than 
one part in 10^. Since the eccentricity of an orbit affects how easily 
it can be approximated. Table (3) shows results for a 12^ satellite of 
different eccentricities. Elliptic orbits were used to study the 
relationship between eccentricity and degree needed to fit. More complete 
results are given in Tables III-VII of. [5], and the information 
in those tables may be used in estimating the degree needed to fit a 
particular orbit. The numbers in parentheses are the degrees needed 
to fit real rather than elliptic orbits. It is seen that the results 
for elliptic orbits carry over very well for low to moderate eccentricities. 
For higher eccentricities, the degrees calculated for elliptic orbits 
are too small, but they still give an idea of the degree needed for 
fitting a real orbit. 

Conclusions and Recoinrendations 

Experience has shown that polynomial approximations in terms of 
Chebyshev polynomials are very effective in representing, satellite and 
planetary ephemerides. They meet all the criteria we specified for 
compression, accuracy, and spanning several orbits, and computationally, 
they are the simplest representation possible. To take full advantage 
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of their capabilities. It Is necessary to choose reference points crowded 
towards the end points of the time Interval. If uniformly spaced reference 
points are used, accuracy Is sacrificed greatly, or an Inordin- 
ately large number of reference points Is required to»find 
that minimum degree. The package will accept conventionally tabulated 
ephemerldes with uniformly space points and Interpolate to obtain the 
proper distinction of reference points. It Is smart enough to abort the 
fitting procedure early when It sees the degree Is too low, and then It 
will automatically Increment the degree and start over. In this mode 
of operation. It Is vastly superior to an L.S. formulation. Even In 
the case where the desired degree Is known, it Is superior to L.S. fitting 
in that the algorithms are very efficient and they automatically give 
the error In the fit. This package Is available In card form from 
the Space Systems Division of NRL. Supplied with It Is a set of test 
data with results which can be used as a benchmark. At present. It 
has been successfully run on Amdahl computers at the University of 
Cincinnati and the Draper Lab, on the NRL Advanced Scientific Computer, 
on the NRL DEC System 10, and on the NBS CDC computer. 
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APPLICATION OF SEMIANALYTICAL 
SATELLITE THEORIES TO PRECISION 
ORBIT DETERMINATION 

P. J, Cefola 

C. S, Draper Laboratory 


Introduction 


Over the last several years, various space mission cen- 
ters including NASA/GSFC, NASA/MSPC, JPL, NORAD/ADCOM, SAMSO/ 
Aerospace and AFGL have supported the develooment of semi- 
analytical satellite theories based on the Method of Averages. 

The intent of all these efforts has been to produce 'fast' 
orbit computational algorithms for mission analysis, mission 
planning, orbit determination and other application programs. 

To date, this effort has concentrated primarily on the de- 
velopment of the equations of motion for the averaged orbital 
elements — that is, the osculating orbital elements minus the 
short-periodic effects. Results include 

— recursive analytical formulations for computing the 
averaged element rates due to gravitational perturba- 
tions (zonals, tesseral-resonance, and lunar-solar 
effects) (see References 1-6)*. 

— the development and refinement of numerical averaging 
concepts for computing the rates of the averaged 
elements (primarily for atmospheric drag) (see References 
8 - 12 ). 

— the widespread application of nonsingular orbital ele- 
ments in formulating the averaged equations of motion 
(References 2,3,5,6,8,16, et. al.). 


• The pioneering wor)c of Hunt Small (R^'ference 7) should bo 
noted in any discussion of recursive satellite theories. 
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— investigations into the efficiency of various numerical 
integration processes for the solution of the averaged 
differential equations of motion (References 9 and 11) . 

— > investigations into various methods for computing the 
averaged elements at a particular epoch, g:<ven a high 
precision state vector at that same epoch (the 'osculating 
to mean' transformation) (References 8, 9, and 11). 

These results have been sufficient to ma)ce *"he semianalytical 
theory the preferred choice when very long data arcs are in- 
volved and when modelling of the short-periodic oscillations 
^s not required. Thus semianalytical theories are usea for 
many of the c rbital computations in preflight mission analysis 
(for example, see Reference 13) . The seune situation holds 
when long arcs of averaged element data are being processed to 
construct geopotential fields or atmospheric density .odels 
(Reference 14) . 

In addition, preliminary consideration has been given to 
the computation of the short-periodic oscillations at the out- 
put points given the averaged elements, in the context of a 
semianalytical theory. Lutsky and Uphoff (Reference 10) pro- 
vided an approach for computing first order short-period ics 
that could be attached to their numerical averaging program. 

Very promising numerical results, with respect to the accuracy, 
are provided in Reference 10. Vashkovjak (Reference 15) provided 
a detailed treatment of the short-periodics for the 24-hour 
synchronous equatorial orbit in tile context of a semianalytical 
theory. Again, high accuracy wa» obtained. And, of course, 
the first transformation of canonicol satellite cheory pro- 
vides the formulas for the recovery of the short-periodics 
due to J2 (Reference 16). 

Despite these efforts, the semianalytical theory has not 
been accepted as a replacement for the Cowell method of Special 
Perturbations in applications where high accuracy output is 
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required frequently (for example, see Reference 17). This 
requirement for frequent output corresponds exactly with the 
requirements of precision orbit determination. 

The outline of the remainder of the paper is as follows. 
First, those factors which limit the usefulness of current 
implementations of the semianalytical approach are discussed. 
Next, numerical and analytical enhancements to the semianalytical 
approach are discussed. Finally, a simple mathematical model 
is provided to estimate the computational speed of a semiana^ 
lytical theory .mploying the suggested enhancements. The model 
can factor in current experience with semianalytical theories 
(integration stepsizes, quadrature orders, speed of recursive 
foirmulaticns, etc.) and the chuacteristics of the particular 
output requirement [observation span (or orbit determination 
interval), observation rate, obsejrvation model, etc.]. Compari- 
sons with numerical integration (Special Perturbations) are 
suggested. 

Evalution of Current Semianalytical Theories vs. the Precision 
Requirement 

To the author's knowledge, there have been only t«ro ser- 
ious studies of a semianalytical theory in a precision orbit 
computation application. These are: (1) the evaluation of 

the KAESTRO numerical averaging theory for a detailed mission 
planning program (Reference 17) and (2) the evaluation of the 
GTDS numerical averaging theory for an crbit determination 
application (Reference 18) . 

Reference 18 concluded: 

— that it was possible to successfully fit the averaged 
dynamics directly to raw observation data 

— that observation editing criteria designed for Special 
Perturbations DC's might lead to the rejection of good 
data in an Averaged DC since short-periodics were not 
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modelled . Even if the editing criteria are relaxed, 
neglect of short-periodics might cause the data to 
appear biased over short observation spans (this is 
because the short-periodic oscillations are much 
larger than typical observation errors) 

— that the multistep numerical integrator did not exhibit 
the full advantage of the semianalytical theory for 

1 or 2 day orbit determination intervals. 

Reference 17 's prime concern was with the computation 
time characteristics relative to Special Perturbations. There 
are two points that seem important to mention with respect 
to this study: 

— that output was required every 2 minutes over a 7200 
minute span (5 days). This requirement was imposed 
in order to simplify (in terms of software changes) 
the interface between the semianalytical theory and 
the application program 

— that the stepsizes employed in the numerical integra- 
tion of the averaged equations of motion were severe- 
ly constrained first, by the retention of the tesseral 
m-daily effects* ** and second, by the use of numerical 
averaging. 

Desirable Enhancements to Current Semia laly tical Theories 

We first list desirable enhaucentsnts to the semianalytical 
satellite. The theory implemented in the R&D version of GTDS 
is ta)cen as the baseline. The enhancements are: 

1) a self-starting low-order integrator with a matching 
Interpolator 

* The 'm dailies’ are due to the terms in the tesseral harmonic 
potential which only depend on the slowly varying satellite 

orbital elements and the Greenwich Hour Angle. 
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2) a recursive analytic formulation of the short-periodics 
Including 

— zonals 

» nonresonant tesserals (including m-dallies) 

— lunar-solar 

3) a low-order interpolation of the approximate high pre- 
cision position and velocity within an observation pass 

The self-starting low-order integrator is intended to 
take advantage of the fact that the integration stepsize for the 
averaged dynamics (1 to 4 days) is large relative to the ob- 
servation spans typically associated with high precision batch 
differential corrections (for example, 1 or 2 days for Landsat) . 
The matching interpolator provides the averaged elements at 
any observation time within the step without accessing the 
averaged force models. A low-order Hermite procedure (References 
19 and 20) may be appropriate. 

The need for an accurate model of short-period ics seems 
obvious in a production orbit determination environment. It 
is fortunate that first order models of 'the short-periodics 
are thought to provide accur icy down to 10m (see Kozai, Ref- 
erence 21) . M-dailies are iacluded in the output time compu- 
tations so as not to constrain the stepsize of the averaged 
integration. 

Since observation rates are on the order of 6 or 10 ob- 
servations per minute and since the grid interval for short- 
periodic interpolation is in the range of 2 to 10 minutes 
(Reference 20), the computation of the ephemeris*at each ob- 
servation time via an interpolation procedure seems to make 
good sense. Thus we will utilize the analytical model of 
8hort-periodic5 only on the interpolation grid and not on 
the much more 'dense' observation time grid. 
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Satellite Theory Simulation Model 


Straightforward analysis of the semianalytical theory 
described in the previous paragraphs leads to the following 
model of the CPU time. 

CPU TIME = (Nmj^ + 1 ) (t2 + ffl3t3) + m^m^tj ( 1 ) 


» number of integration steps 

N >> function evaluations per step 

t2 “ time for analytic contribution to single-averaged 
element rates 

m3 «= number of density evaluations per quadrature 

t3 * time for each evaluation of the quadrature integrand 
* number of observation passes 

m^ * number of output points per pass which require analytic 
short-periodics 

t^ »* time for each computation of .the analytic short- 
periodics 

Note that Eq. ( 1 ) concentrates on the high cost mathema- 
tical functions. Eq. ( 1 ) does not attempt to model the various 
interpolation procedures although it does include the generation 
of the data required to constr^^^the interpolation coefficients. 

Table I provides sample evaluations of Eq. ( 1 ) for two 
typical scenarios. In both cases, it appears that significant 
computational advantage can be obtained via the semianalytical 
method. This is because we expect t2 and t^ to be on the order 
of a high precision perturbing acceleration evaluation. This 
has been demonstrated for t2 in Reference 4 . For t^, this 
expectation is based on mathematical analysis performed to date. 

However, this advantage is dependent on the enhancements 
described in the previous sections. For example, suppose in 
Case 1 that the m-daily effects were retained in the averaged 
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Table I - Typical Scenarios 


Case 1 - Very Low Altitude 

1 day arc 

16 rev/day orbit 

2 lainuteP/'obs pass 
6 passes/day 

3.0 cbs/minute 

= 1 step 
N - 4 
nij « 24 

«4 * 6 

ittg = 2 

CPU « 5 tj + 120t3 + 12tg 


Case 2 - Medium Altitude 

2 day arc 
14 rev/day orbit 

5 minutes/obs pasp 

6 passes/day 
10 obs/minute 

nij^ * 1 step 
N = 4 
Blj * 9 

n»4 * 12 

®5 “ 3 

CPU - St2 + 45t3 + 36tg 


REMARKS: 


RFMARKS; 


1. Cowell orbit generation would 
require around 1600 steps. 

2. Total obs “ 120/arc 


1. Cowell orbit genera- 
tion would require 
around 2800 steps. 

2. Total obs = 600/arc 


integration. Then might grow to 8 or 10 steps for the 1 
day arc. For = 10, Eq. (1) gives (for Case 1) 

CPU = 41t2 + 984t3 + 12tg (2) 

Clearly most of the advantage would be lost for all reasonable 
models of the atmospheric density**. Or suppose that the multi- 
step numerical integration procedure was retained. The starter 


** This corresponds with the configuration of the semianalytical 
described in Reference 17. 



associated with this process would reduce the advantage of 
the seinianalytical approach. In Case 2 , the advantage clearly 
depends on offloading the computation of short-period ics from 
the grid of observation times to the interpolation grid 
(3 points/pass) . 


Conclusion 


While the above simulation exercise suggests that the 
aemianalytical approach can be very desirable, what is really 
needed is the development of an OD test-bed employing this 
approach. Such a test-bed could demonstrate the advantages 
and disadvantages in an unequivocal manner against real 
observation data and scenarios. 
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CRITICAL INCLINATIONS IN SATELLITE THEORY 
Andre Deprlt* 

Department of Mathematical Sciences 
University of Cincinnati 

The main problem of satellite theory is described in polar 
coordinates by the Hamiltonian function 

^ " 4,0 " ■‘^2’ 


T/o = 


^'- 2 ) - ^ 

r r 


n/ _ u 3 2, 3 2 - T 

” r 2 ^ 2 " 4^ " 4^ 20 ]. 

It is proposed to find a solution of it with the following 
properties: 

1°) the reference orbit is Keplerian; 

2°) no restriction is imposed on the eccentricity; in 
particular, it is exempt of singularities - real or 
apparent - for small eccentricities; 

3°) no restriction is imposed on the inclination; in parti- 
cular, it is exempt of singularities - real or apparent 
for small inclinations; also it is valid even in the 
neighborhood of inclinations at which the perigee is 
stationary. 


On leave at the Division of Applied Mathematics, National Bureau 
of Standards, Washington, D.C., 20234. 
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The construction proceeds in two steps. 


In the first step, a can^ical mapping, called the elimination of 
the (equatorial ) parallax , changes the Hamiltonian into the function 



, r 1 n /OI\ 

^ = I ~ t (-) 


a,2n 


n>0 n! 


4- = y I 

^ OiiiLn/2] 0<j< 


I 

0<k<n 




9 


X * e COS g, 
Y = e sin g, 
s = sin I. 


In the second step, a canonical mapping, called the revolution 
in orbit , changes the Hamiltonian into that of a Keplerian system. 

Both transformations are obtained in application of a perturbation 
algorithm based on Lie Series. The basic differential equation 


lltV) 

^^0’ n^ ^0 


may be reduced to an elementary quadrature if one makes the following 
observations. 

7.,-> 

i) Assume that</ is of the form 
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- 


Then 






where designates the partial derivative of with respect to 
its first arrument (namely e). 

(1i) Throughout the construction of the Lie triangle, the elements 
be maintained in the form 

J 

J r p J 

Therefore, at the end of any row in the Lie triangle, the partial 
differential equation reduces to the quadrature 

»i“„ = «S - "S- 

In the course of eliminatirq the parallax, the factor Hq emerge as 

finite Fourier sums in the argument 0 of the latitude. It is thus 

natural to set the unknown factor Hq equal to the average of ]ltQ. Hence 

will be a purely periodic function of 6. 

In the second transformation, the unknown factors Hq are set to 

zero. Hence will be a finite sum of mixed terms 0^ sin je and 
n 

0^ COS J0. 

At the first order in e, the revolution in orbit transforms the 
argument of latitude according to equation 
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e *- e'[l + - 5c'^)]. 


So the rotation of the coordinate system implied by the canonical 

mapping becomes the identity at the inclination of stationary 

-1 2 

perigees, namely I * tan 2 for which 1 - 5c =0. This explains 
why the solution does not recognize the inclination of stationary 
perigees as a critical singularity: the revolution in orbit adjusts 

the frame of reference so that it follows the perigee. The property 
is typical of a non-essential resonance of type (1:1) whereby a 
rotation of th<» coordinate axes preempts the apparition of zero 
divisors. 

The calculations have been executed by hand - with the collaboration 
of Mrs. Deprit-Bartholome - to order 2 for the elimination of the 
parallax and^^ order 3 for the revolution in orbit. The results 
have served to check uie computer programs which then carried out 
both transformations to order 4. 



Th^new theory is the first one to have obtained the fourth order 
rms. The most accurate observations currently available require 
that the main problem of artificial satellite be solved to order 3. 
The terms of order 4 will serve to estimate the errors induced by 

3 

truncating the series beyond t . 

The generating functions for both mappings are much smaller by 
the number of terms than those of the conventional (Kozai) and not 
so conventional (Aksnes) theories. 
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The computer programs to execute the reduction in a literal 
form involve a processor for Poisson series. The latter is the 
latest version of MAO. From a package of subroutines written in 
assanbler or in Fortran, MAO has evolved into a subroutine gen- 
erator. At compilation, macro variables are set up to specify 
the type of Poisson series needed to solve a particular class of 
problems. The generator is coded to be preprocessed and compiled 
by the IBM optimizing compiler for PL/I. It will be made available 
to mathematicians In dynamical astronomy and non-linear mechanics 
as soon as the documentation has been published. The Department 
of Astrontxny at the University of Thessaloniki is considering 
transferring MAO from IBM to UNIVAC. 

In the course of expanding the functions generating the 
canonical mappings to solve the main problem of satellite theory, 
a "profiler" in line counted how many times the subroutines in the 
package were called. There have been 9452 algebraic and differential 
operations on Poisson series, 286773 "list" operations (to find or 
to create nodes in chains) and 216651 algebraic operations involving 
rational numbers (represented and maintained as pairs of relatively 
prime decimal integers). The execution time was 70". 54 on an 
Amdahl 470-4 operating under OS/VS-2 at the University of Cincinnati. 

St>uuye 6poAX! WheAe. dutincution has no place. 

Ok none, and may be anyukcAe loe choose! 

WheAe MAO, comlXted to kU endless nace. 

Runs tike a madman diving {pA iXs Aepose! 
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A SINGULARITY FREE ANALYTICAL SOLUTION OP ARTIFICIAL 
SATELLITE MOTION WITH DRAG 

Alan Mueller 

Analytical and Computational Mathematics, Inc. 
1275 Space Park Drive, Suite 114, 

Houston , Texas 77058 


If an analytical satellite theory which includes the drag 
perturbation is to be successful, it must have three important 
qualities. First the theory should be based on a canonical 
formulism whereby one can use the powerful tools provided by 
hamiltonian mechanics. Secondly, the model used to describe 
the forces acting on the satellite must not be so simplified 
that the theory becomes only a mathematical exercise. Lastly, 
the resulting theory must be concise so that the accuracy 
gained outweighs the extra computer costs required to reach 
that accuracy. 

Scheifele (reference 1) has developed an analytical sat- 
ellite theory based on the regular, canonical Poincare-Similar 
(PS<|>) elements. This is a very powerful set of elements which 
are in an extended phase space and have an independent variable 
which is similar to the true anomaly instead of time (references 
2, 3 and 4). A very accurate and concise satellite theory 
has been developed to include the first order short period 
and secular perturbations of an oblate central body. The drag 
theory has been built on top of the J 2 theory. 

The assumption in this theory is that the drag force is 
tangential to the orbit and proportional to the square of the 
velocity magnitude of the spacecraft. The constant of pro- 
portionality, which is a product of the density of the atmo- 
sphere, the ballistic number, and the drag coefficient, was 
not specified in Scheifele's theory . Since the lifting force 
relative to the drag force and the inertial velocity of the 
atmosphere relative to the satellite velocity is small, the 



model used is adequate for giving the direction of the retard- 
ing force due to the atmosphere. Thus an important contribu- 
tion to the analytical solution was made. The report (refer- 
ence 1) is a concentrated e'ffort to canonically transform the 
forces into the PS space and also place them in a form suit- 
able for solution. Therefore, the direction of the PS canon- 
ical forces has been determined but their magnitude was not 
completely specified. Also, the tools of hamiltonian mechan- 
ics were used to transform the forces correctly and reduce 
the size of the equations. Due to the character of the PS 
system, the equations which describe the motion are relatively 
simple and thus the first and third qualities mentioned above 
are satisfied. 

The second task was to develop an atmospheric density 
model that can be used in Scheifele's theory. In developing 
a density model for the analytical theory one is severely 
restricted by the fact that the model must be in the form of 
a fourier series in the true longitude. As is the case in 
most analytical theories, the perturbations must be written 
in a fourier series to facilitate solution by quadrature. 
Several density models have been developed to predict very 
accurately the density at any point in space and time. Ex- 
amples are the Jacchia model (reference 5) and the USSR model 
(reference 6). Both models are extremely complicated and too 
unwieldy for analytical satellite theories. In the analytical 
theory of Brouwer and Hori (reference 7) the density model was 
assumed to be an exponential function of the radius. However, 
the atmopsheric density is strongly dependent on the sun and 
its position, and also the oblate figure of the earth. Thus 
Brouwer's assumption is not valid. A completely new model 
needed to be developed which is both accurate and can be 
written in a fourier series. 

In developing the new model, the approach taken was to 
construct a model which is able to simulate the Jacchia density 
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along a particular orbit. The value of the coefficients in 
the new model are determined by a procedure called **calibra- 
tion". A simple formulation allows the model to be inverted, 
i.e. given the density at different points along the orbit 
(as determined from Jacchia), one can compute the coefficients 
of the model. The coefficients are implicit functions only 
of long period effects and can be conjidered constants in the 
analytical theory . 

The model has been fit to a particular orbit to include 
the variations in the observed density due to two-body changes 
in the height, and the two-body changes of the angle between 
the sun and the satellite (diurnal effect). Included in a 
manner similar to that of Santora (reference 8), are the 
density variations caused by changes ii height due to the ob- 
late figure of the earth and the snort periodic oscillations 
in the radius due to . The density model also ’’accounts*’ 
for the changes in the density because of secular perturba- 
tions in the height due to drag itself. 

The result is an accurate density model which can be 
implemented into the drag theory. Numerical experiments de- 
monstrate the close agreement between the new model and the 
Jacchia model. 


The last stage of the analytical theory is under develop- 
ment. This involves constructing a computationally efficient 
manner J^^uich to expand the equations of motion into a 
fourier^lCTies . This requires a careful balance of explicit 
manual computation, explicit equat^P^s by computer manipula- 
tions, and lastly but not least, the recursive relations. 


Most, but not all of the theory has been implemented on 
the computer. Comparisons with numerically integrated solu- 
tions verify that the analytical theory is extremely accurate. 
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Effect of Atmosphere on Venus Orbiter Navigation 


Mohan Ananda and Braulio Sanchez 


Jet Propulsion Laborator5' 
Pasadena, California 



ABSTRA CT 

The current uncertainty for atmospheric models of Venus Is significantly 
large. The orbital prediction requirements for Pioneer Venus Orbiter with 
Its relatively low perlapsis altitude (150 km) have brought concern on 
navigation capabilities. This paper investigates simplified but realistic models 
of the Venusian atmosphere on orbit determination accuracy. A model with 
polynomial representation of the atmospheric scale heights Is assumed for 
statistical error analysis. Covariance analyses have shown the effect of model 
errors in the Venusian atmosphere can be minimized for trajectory prediction 
after processing several orbits of data. Studies include the sensitivity of 
periapsls data, arc length, data rate and station coverage for determining 
atmospheric parameters. Periapsls data are highly sensitive to the gravity 
field. The gravity field of Venus is essentially unknown and thus it is 
necessary to determine both gravity and atmospheric parameters simultaneously. 
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DENSITY MODELS TOR THE UPPER 
Douglas L. Dowd^ and B. B. 


atmosphere'^ 

Tap ley ^ 


1 . Introduction 

One of the more important problems associated with the task of defin- 
ing the orbit of a near earth satellite is that of modeling t.he effects of 
atmospheric drag. Errors in the drag model can lead to significant errors in 
the determination and prediction of the satellite position. The drag acceler- 
ation is modeled by Che relation 






V 

r 


where p is the atmospheric density, is the drag coefficient, A is 

the cross sectional area normal to the relative velocity vector, m is the 
satellite mass and is the velocity vector relative to the atmosphere. 

Hence, the uncertainty in the drag acceleration can be separated into three 
components: a) errors in the atmospheric density model, b) errors in the 

ballistic coefficients, and c) errors in the satellite relative velocity- 
The first of these error sources is due to inaccuracies in a priori models 
and presents a limiting factor in the accuracy with which the velocity and 
position of an orbiting satellite can be determined. 


^Aerospace Engineer, Mission Planning and Analysis Division, L. B. Johnson 
Space Center, Houston, Texas. 

2 

Professor, Department of Aerospace Engineering and Engineering Mechanics, 
The University of Texas at Austin, Austin, Texas. 

X 
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Xorn^ally, the atmospheric density is modeled by defining; an a priori 
static model based on historical satellite tracking data. Since the atmo- 
spheric density depends on such external influences as solar and geomagnetic 
activity, computed values of the density will be in error due to inaccuracies 
in the original uefinition of the density model as well as time lags in updat- 
ing the parameters which account for the effects of solar and geomagnetic 
activity. 

In a number of contemporary satellite missions, the requirement for 
performing the orbit determination and prediction in real-time has placed 
an emphasis on models which, in addition to being accurate, require a 
minimum of computation time. In addition, if the computations are to be 
performed using a satellite-borae computer, the models must be efficient with 
re;'.arJ to computer storage requirements. 

In this investigation, consideration is given to three contompoiary 
atmospheric density models which have been selected as the test candidates 
to meet these requirements. The models considered are the Analytic Jacchia- 
Koherts Mode] [1], the MOiliflcvI Harr is-Pr ios cor Model [2], and the U.S.S.K. 
Cosmos Satellite Derived Density Model, corarnonly known as either the Russian 
Model or the U.S.S.R. Model [3]. Each of the models and their respective 
variacions is discussed separately, and a comparison of the computational 
characteristics of Che models is presentee. Finally, recommended modifica- 
tions for improving both the computation speed and accuracy are presented. 

2. The Analytic Jacchia-Roberts Model 

The Analytic Jacchia-Roborts Model calculates atmospheric densities 
for altitudes at 90 km and above. The model, an analytic representation of 
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1970 tabular density r«:odal [4} developed by Roberta (1] incorpotatta 
tiia revisions to the tubular tnixicX which were published in 1971 by Jacchia {5J. 
The oethod divides the upper atmosphere into three altitude bunds for the 
calculation of atmospheric density* Specifically, these bands are 90-100 km, 
100-125 km, and higher than 125 km. The terminal conditions in each lover 
band are the initital conditions for the nesit higher band. ’*'hercfore, the 
dctcnnlnnCion of the density within any j»*von altitude baud roqulrcs the 
calculation ^of the terminal conditions in each of the lover bands. The method 
is predicated on the assumption of a termperature profile and assumed values 
for the molecular mass of the major atmospheric constituents. The atmospheric 
density Is dctcmiincd then by integrating cither the barometric equ.'itioa for 
altitudes from 90 km to 100 km or the diffusion equation for altitudes above 
100 km. The major constituents considered by Jacchia are nitrogen (^ 2 ), 
argon (Ar), helium (He), molecular oxygen ( 02 ), atomic oxygen (0), and 
hydrogen (H) . 

For altitudes in excess of 125 km, the termperature profile is de- 
fined .-m; thematically by an asymptotic function. Jacchia originally chose to 
use the inverse tangent function vhich did not produce an exact differential 
In ihe diffusion equation and hLs tabular model is a result of numerical 
integration of t)ie diffusion equation. Roberts [1) assumed an exponential 
temperature profile which allows for the analytic integration of the diffusion 
op*-*ration. In oUlier of the assumed topiper.iture profiles, there is no mathe- 
matical upper altitude limit. As altitude increases without bound, the 
te.T.pcrature asymptotically approaches the exospheric temperature. At some 
unspecified altitude, the density of the atmosphere has decreased to the 
point where the gas atoms move in ballistic trajectories and no longer 
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Interact with one another in support of the laws of fluid dynamics. There** 
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fore, temperature profile approaches its bound, the validity of the 

diffusior. becomes suspect. The altitudes at which this occurs in the model 
is highly dopeudent upon the value of tb:^ exospheric temperature [5] and 
ranges from S80 km for an exospheric temperature of SOO^K to 2000 km at 
1900‘^K. 


In the original Jacchia model, and in the analytic model as well, the 
exospheric temperature (T^) calculation accounts for the observed variations 
in the density. The variations in T^ are correlated with variations in the 
90-day mean flux of the solar 10.7 cm radiation, where 


10.7 


Ci 


rT 

i F. 


10.7 


dt ] X 10 


-22 -2 „ -1 
m Hz 


and with the daily variations of y from the mean. The value of the 

y solar flux varies with an 11-year period while the F^^ y flux has 
a 27-day period with an apparently random amplitude due to the effects of the 
solar rotation. The temperature calculation also accounts for the variation 
in density as a function of the local solar time of the point in question 
^'di’urr.al effect) and changes in the geomagnetic activity. The atmospheric 
density determined from this exospheric temperature is corrected for the 
seasonal-latitude variations of helium, and the variations in hydrogen 
concentrations above 500 km. A simple logic flow chart of the Analytic 
Jicchia-Robcrts Model is shown in Figure 1 while the specific algorithm for 
the atmospheric density model is described in detail in Reference [6]. 

An efficient modification to the basic Jacchia-Roberts Model has 
b».en adopted for use in the Goddard Trajectory Determination Subsystem 



(CTijS) 12], The cs;sential modification^ are that the atmospheric density 
at ICO km. and the density numbers of the ro«ijor atmospheric constituents at 
125 km. arc all approximated by sixth deijree polynomials in Instead of 
being calculated as the terminal conditions of the two lowest altitude 
bands. Further discussion regarding the computational aspects of these 
modifications is given in Section 6. 

3 . The Modified Hnrrls-Prios tor Model 

The Modified Harris-Pricstcr Model [2] is based on an extensive 
tabular utaiic model of the uiipor atmosphere in the altitude band from 120 
km. to 800 km. [7]. The first modification of the Harris-Pricstcr Model, 
accomplished at the NASA Goddard Space Flight Center, was to exponentially 
extrapolate the tables to include altitudes down to 100 km and up to 1,000 km. 
The model, as incorporated into GTDS [2], retains its tabular form in a modified 
format. In the original formulation, there are 10 separate tables, each being 
associated with a particular value of the smoothed (5-month average) flux of 
the solar 10.7 cm radiation r.. aging from ^ 65 to ^ • 275. This 

ranp,c of ^ encompai»bes tlio total variation of F^^ ^ over the 11-year 

solar cycle. Each table consists of 12 subtables which list the atmospheric 
densities for the local solar time at 2-hi»ur intervals. The tables u>r tlie 
Modified ilarris-Priestcr Mcd^l are formetl by extracting the maximum and minimum 
densities for each altitude fi\)«i the bub tables for each' solar flux level. The 
absolute maximum and minimum are chosen without regard for the local solar 
time. This was done because the diurnal maximum and minimum densities do not 
a 'pear in the tables at 1400 hours and 0400 hours at altitudes below 320 km, 
as is the case for the observed extrema. The Modified Harr is-Pries ter Model 


42 



then derives the straosphLrlc density fron a set of 10 tables associated with 
the smoothed flux of the 10.7 cm *{olar radiation , where each tabic relates a 
diurnal maxJmina nnd minimum density for each of the tabulated altitudes from 
100 to 1,000 kiu. 

The atmospheric density for a given altitude is determined then by 
entering the table as Delated with the value of F ^ - most nearly equal to 
the measured solar flux, exponentially interpolating the maximum and minimum 
densities with respect to altitude, and then applying a cosine interpolation 
for the diurnal variation. This procedure yields a density distribution which 
is symmetric with respect to the apex of the diurnal bulge. The apex of 
the diurnal bulge Is assumed to follow tlu.* subsolar point by 30® in the same 

latitude. It Js known (8] thac li*e observed diurnal variation is not 
sy4.mietric and tiic Analytic Jncchia-Robcrts [1] and the U.S.S.R. Cosmos Satel- 
lite Derived Models [3] account for the asymmetry. The Jacchla-Roborts Model 
iiccomplishea this by computing an asymmetric tomporaturo distribution fr%Dra 
which density is determined. In this investigation, a similar procedure has 
been applied directly to the density computation to provide an asymmetric, 
density distribution in the Modified Harris-Prlester Model. The model with 
this procedure is referred to as the Asyrametric Modified Harris-Prlester 
Model. The detailed computational algorithms arc given in Reference [6]. 

The pertinent equations for describing the diurnal variation in the 
Modified and As^nnmctrlc Modified Harris-Prlester Models are as follows: 

Modified Harris-Prlester Model 
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r 

% 


wh«r« h la tha altituda. and p^ ^^ (h) ara tha intarpolatad 

daily Blnlmua and naxlmum danaltlaa from tha imdlflad Cablaa [2], and ip 
la tha angla batwaan tha gaocuntrlc poaltion vaetora of tha point whara tha 
nodalad danaity ia daairad and tha apax of tha diurnal bulga< 


Asymmetric Modified Harrls~Prieater Model 

p(h) - Pjj(h) + [ Pjj(h) - Pjj(h) ] coa“ (x/2) (3) 


where 


Oj,(h) 

Pp(h) 

0 

n 

T 


- I !♦ + «l 

- y !♦ - fi| 

■ H + 6 + A sin (H + y) ("if < t < ir) . 


(4) 


In the above equations, H is the local solar time, d la the 
geographic latitude of the subsatelllte point, and d ia the aolar declln- 

4 w • . . Vt..4.UCw »*iC ' k. 1. 4- w XI* t. XOn,:? 

from Jacchla'g temperature equation [3 ], are: 
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The shape of the diurnal bulge » as modeled by the Asyometrlc Modified 

Harris^Priester Model, is illustrated in the polar plots in Figures 2 and 3* 

In Figure 2, the angle measure is latitude, and the magnitude in the radial 

deviation represents the noimalizcd modeled density variation at some assumed 

constant altitude h where p , (h) » 1.2 and p (h) • 2*0 are the 

mxn max 

assumed der4sity values at h. The specified solar hour angles H are for 
Che right halves of the plots with the hour angle for the left halves being 
H + 180^. In both Figures 2 and 3, the assumed solar declination is 15*. 

The unit circles in each figure are included to emphasize the changes in 
the density magnitude. In Figure 3, the angle measure is longitude 
(or solar hour angle) measured from noon, and the radius magnitude repre- 
sents the modeled variation at constant altitude and latitude. The 
Figures 2 and 3 show the global maximum density occurring at the subsolar 
latitude and 31,226'* after noon and the global minimum occurring at 

latitude and 137.01® before noon. 

0 

A . The U.S. S .R. Cosmos Satellite Derived Density Model 

Ine U.S. S.R. Cosmos Satellite Derived Density Model is based on the 
observations of 145 Cosmos satellite«i over the time period from 1964 through 
1970 [4]. The model determines the atmosp/.eric density directly by substi- 
tuting the input parameters into a set of equations containing twenty coeffi- 
cients derived by fitting density observations over the range of altitudes 
and temperatures encountered by the ICosmos satellites. The use of the 
current model is restricted because the coefficients were empirically deter- 
mined over a limited altitude region and during only a portion of the 11-year 
solar cycle. The data were extended by using Jacchia's 1970 Model, but the 
altitudes for which the ir.od< * valid is still only between 140 and 500 km. 
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I The coexilcients are given in four sets for four reference values of the 

I «22 

I 10.7 cm, solar radiation flux; specifically, 75, 100, 125, and 150 x 10 

I 

I —1 

\ WATTS m Hz • Since the model uses the reference value which Is near^ 

j esc the 6-montb average of the daily F^^ ^ , the model is valid for 6<-month 

j averages of F^q ^ between 65 and 165 x 10^^^ WATTS m ^ Hz 

I . The details of the U.S.S.R. Model are given in Reference [6]. The 

I essence of the model is chat the nighttime density profile is corrected 

I 

^ by four multiplicative factors K^, i • 1,2»3,4. The K^- factors 

I include corrections for the diurnal density variation, , the dally 

J variation of F^^ ^ , the observed semi-annual variation in density, 

j , and fluctuations in geoMe;aetic activity, . 

The total density is ihc,\ represented by the equation 

r 

J p “ ph h h h 

5. Explanation of Atmospheric Density Profiles 

A study of atmospheric density and its effect on the motion of a 
near earth satellite would not be complete without an explanation of the 
correlation between the orbit of the satellite and the density profile which Is 
encountered by the satellite. Normally, the atmosphere is discussed as a 
separate system with density presented as a function of altitude for various 
values of the other parameters vULch have been correlated to variations in 
the observed densities. In this discussion, the intersection and interaction 
of two dynamical systems, the atmosphere and the orbiting satellite, will be 
considered. The density profiles which will be discussed are referred to 
f the Modified Analytic Jacchia-Roberts Model since this model contains all 

of the essential variations wliilo retaining computational ofricloncy. 
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It is not difficult to understand the relationship between atmospher- 
ic density and altitude— as altitude above the earth* s surface increases t 
density decreases, provided everything else is constant. Therefore* if the 
variation of the satellite altitude is known as it moves in its orbit, one 
would expect to see an inverse variation in atmospheric density. Obviously, 
orbital eccentricity has a large effect on the altitude variation. Con- 
sidering that the earth is not spherical, orbital inclination also has an 
influence on the altitude variation as does the orbital perturbations due 
to the uonsphericity of the geopotential. To illustrate these points, 
refer to Figure 4 which shows time histories of altitude above the refer- 
ence ellipsoid, geocentric radius and atmospheric density for one orbital 
period. The orbit used to /.oiUTvile these results is approximately circular 
with initial osculating Keplerian elements as follows 

a « 6682473 m 
e « .000646254 

i - 67. 99® 

w - 0.0*^ 

a - 0 . 0 * 

E - 0.0® 

16 *^ 2 '' IVC. ]<)73 

4 .C should be noticed that the amplitudes of the variations in radius and 
altitude arc not of the same magnitude which indicates the dual effect of 
the earth's nonsphericity on the altitude variation and, in turn, on the 
density profile. The density curve indicates chat there are other major 
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effects involved In shaping the density profile. To aid in the identifica- 
tion of the Jtosc obvious of these effects* consider the histograms in 
Figure 5 which were generated exactly as those in Figure 4* except that the 
nodal line has been rotated ninety degrees* i.e. H « 90*. The altitude varia- 
tion and latitudinal displacenKsnt from the diurnal bulge are the same in both 
cases* whereas the longitudinal displacement from the diurnal bulge is offset 
by ninety degrees. There is a marked difference between the density profile 
which appears as a phase shift in an apparent once-per-revolution variation. 
This differerence illustrates the diurnal variation and its relative im- 
portance in modeling atmospheric density. 

Up to this point the discussion has related to the shape of the 
density profile. The magnitude of the atmospheric density exhibits other 

I 

variations which are still present when altitude and diurnal variations are 
eliminated. Most significant are the variations in density due to varia- 
tions in solar radiation and the interact Ion of the solar wind with the 
earth’s magnetic field. Density profiles are In Figures 6 through 

8 which illustrate the changes in density that are correlated with both long 
and short term variations in solar extreme ultraviolet (EUV) radiation as 
evidenced by the flux of the 10.7 cm. solar radiation. The effects of 
geomagnetic heating on the magnitude of density are shown in Figure 9 in 
which density profiles arc presented for four values of the planetary 

geomagnetic index K • The initial conditions used to generate the orbits 
P 

for Figures 6 through 9 arc the same as those used for Figure 4. 

6. Comparison of the Density Models 

Each of the density models discussed in the preceding sections will 
yield a density profile along nny given trajectory which is different Chan 


48 



the density computed by any other model. Comparisons of the density profiles 
are shown plctorlally In curves of density versus time in this section. The 
trajectories were generated from the initial osculating elements: 


a " 

6678155 m. 

U) • 

3.02* 

e ■ 

0.0 

n • 

254.26 

i - 

67.99* 

E • 

356.98 


epoch - 16** 2** 47® 55.537® Dec. 1973 


The Nevt^lan equations of motion were numerically Integrated by a fixed 
Step size third order Runge*-Kutta method with a ten second step size. 

A single trajectory was generated and the various atmospheric density 
models were evaluated on this common trajectory. The force model for 
drag used densities from the Modified Analytic Jacchla^Roberts Model In 
the generation of the comparison orbit. The remaining modeled forces used 
to generate the comparison orbit were: 

Two body - u « 3.986013 x 10^^ m^/sec^ 

Nonspherlcal Earth 1969 Smithsonian Standard Earth II to 4th 

Order and 4th Degree 

n body - Solar and lunar gravitational perturbations 

based on Jet Propulsion Laboratory Development 
Ephcmcrls Number 19 


The density profiles shown In Figure 10 are those which would be 
modeled by various forms of the Modified Harrls-Prlester Model. The refer** 
ence profile is the generating density profile modeled by the Modified 
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Jacchla-Robercs Model where 


F 


10.7 


75 


and 


10.7 


75 


K - 1 

P 

The other four profiles shown in Figure 10 are those density profiles which 
were computed by either the Modified Harris-Priester or Asymmetric Modified 
Harris-Priester Models associated with ^ ■ 75 where the shapes of 

the profiles shown were determined by setting the value of n in Equations 

2 and 3 to either 3 or 6. The key to the symbols used to identify the 
curves in Figures 10 and 11 is given in Table 1. 


Table 1. Key to Symbols in Figures 10 and 11 


SYMBOL 

DEFINITION 

MJR A 

Modified Jacchia-Robercs Analytical Model (Reference Model) 

MHP 0 

Modified Uarrifc;-Pricster Model n ■ 3 

MHP - 

Modified Harris-Priester Model n « 6 

^MHP A 

Asymmetric Modified Harris-Priester Model n ■ 3 

AMHP D 

Asymmetric Modified Harris-Priester Model n ■ 6 
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The density profiles shown In Figures 11 are similar to those in 
Figure 10 except that the values of ^ and F^^ ^ used In the models Is 
275 » a value representing the maximum extreme In the ll*year solar cycle* 

Close Inspection reveals that the density profiles determined' by the Asym- 
metric Modified Harrls-Priester Model with n » 3 most closely approximates 
the Jacchla-Roberts profile in shape* It appears that by Judiciously 
scaling and In Equation 4 and by applying a small correction 

to n near the value n « 3 » the Asymmetric Modified Harrls-Priester 
profile could be made tfo very nearly coincide with the Modified Jacchla- 
Roberts profile* 

There Is a systematic difference In the density profiles generated 
by the Jacchia-Roberts and Asymmetric Modified Harrls-Prlester Models 
which Is not apparent in Figure 10* This difference Is a result of the assump- 
tion by Jacchia [4] of a static temperature profile, whereas Harris and Prlcster 
used a dynamic temperature profile to generate their atmospheric density tables 
[7]* These differing approaches are manifested in the models through the 
temperature equations 

Tjj - Tj, (1 + R cos”* n) (2.1) 

Tjj - (1 + R sin”* 0) 

In the Analytic Jacchla-Robcrts Model and the density equations 

• PmIN ^ 

In the Asymmetric Modified Harris-Priester Model. The quantity R appears 
ee a constant in the former model whereas Q is given by 
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W*v, 


which ia not a constant valued quantlty» in the latter model* The effect of 
this systematic difference in the two models is illustrated in Figure 12 
which shows the density profiles generated by the Modified Analytic Jacchia- 
Roberts and Asymmetric Modified Harris^Friester Models during one orbital 
period along a trajectory with initial osculating Keplerian elements 
a » 6662473.58 meters 
e - .00064625 


i * 68.0 degrees 

The difference is most apparent between 3000 and 4000 seconds after the 
beginning of the orbit propagation. In this regard, the Asymctric Modi- 
fied Harr is-Pr tester Model more accurately reflects the real world diurnal 
density variation than the Jacchia-Roberts model. 

Density profiles calcuiated by the U.S.S.R. Cosmos Satellite Derived 
Doaslty .iro coiiip.irid wiili tlic Modified Ana Lytic Jacchia-Koberts prof U os 

in Figures 13 and 14. The initial conditions for the generation of the com- 
p.»ri*;oa <»r!u* t aro iho 1 4«oso u:;od lor Figure 10. The lu'rtinent 

parameters supplied to the models * a generate the prof ilos in 13 




were j » 79.1 , ^ 84.2 , and ^ Jaccra|P«bberts 

Model and ^ « 79.1 , • 75 , and • 4 for the U.S.S.R* model. 

It should be noted tha » 1 and a^ •• 4 are equivalent measures of 

geomagnetic activity. For the profiles shown in Figure 14, the defining 
parameters are F^^ ^ ^ « 150 and « 1 for the Jacchia-Roberts 

profile and ^ • 150 and a^ • 4 for the U.S.S.R. profile. The 

density profiles in these cases are similar in some sense, but not as much 
so as the Asymmetric Modified ii..r/ls-Priester profiles. 





r 
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Comparison of the time required to compute the model densities is 
shown in Table 2 * In terms of computational speed, the U.S.S.R. is the most 
efficient model, with the Modified Harris-Priester Model following closely. 

It should be noted that there is practically no difference in the computation 
time between the Modified and Asymmetric Modified Harris-Priester Models. 

The Analytic Jacchia-Roberts Model requires much more time than either 
the Harris-Priester or U.S.S.K. models. Even though the modification of the 
Jacchia-Roberts Model proposed in Reference 7 reduces the time requirements 
by over 20%, the ocher models arc still more than twice as fast. 

In an independent study, Botbyl [9] investigated the sensitivity of 
the density calculation to the evaluation of a fifth order polynomial occur- 
ring in the Analytic Jacchia-Roberts algorithm. Botbyl shoved that perturbing 
the coefficient of the fifth order term by 1 in the 14th digit resulted in 
density calculations accurate to no more than two or three digits. However, 
the reference density used for comparison did not consider the errors due 
to inaccurate determination of the roots discussed above. To resolve the 
question of the computational accuracy of the model, three density-vs-altitude 
profiles were determined with (1) all single precision arithmetic, (2) double 
precision computation of the polynomial and single precision arithmetic other- 
wise, and (3) all double precision arithmetic. The density digits for the 
three profiles are shown in Table 4. In general, it can be seen that the 


single precision computation is accurate to three or four digits and that 
computing only the polynomial with double precision arithmetic does not 



I 
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•ignificantly increasa the accuracy of the computation. The Importance of 
the increase in accuracy achieved by performing the density calculations in 
dnuble precision must be weighed against the increase in computation time 
and core storage requirements. 

It has been mentioned previously that the Modified Analytic Jacchia- 
Roberts Model is approximately 20Z faster than the unmodified version. A 
comparison of the densities calculated by the* unmodified model using double 
precision arithmetic with densities calculated by the modified model using 
single precision arithmetic is shown in Table 5. The densities were computed 
for T^ from 800^K to 2000^K. The largest error encountered was .082 which 
occurred at an altitude of 125 km. when T^ ■ 2000*K. Between the altitude 
of 90 and 100 km.> the unmodified and modified algorithms are identical. 

The Modified Analytic Jacchia*Roberts Model then is at least as accurate » 
and most of the time is more accurate than the basic algorithm when using 
single precision arithmetic. 

The Modified Harrls-Prlester Model is a very simple, straightforward 
method which displays no computational idiosyncracies. However, potential 
users of the model should consider the physical interpretation of Ls^^ation 3 
which shapes the density profile with respect to the diurnal density varla~ 
ti(‘* * value of the exponent n could bo any real number. Coimnon sense, 

however, tells us that certain values of n would produce profiles that 
almost certainly are not physically realizable. It is not improbable that 
the density variation due to diurnal heating Is a smooth process, that is 
to say at least continuously differentiable. The diurnal variation given 
by the Modified Harris-*Priester Model would then be smooth when n > 1 so 
that n • 1 is an absolute lower bound. However, for values 1 < n ^ 2 . 



the modeled diurnal variation would be such that the profile is broader 
around the maximum than the minimum and this characteristic opposes the 
observed character of the diurnal variation fS]. Conversely* when n is 
large (greater than 8), the density profile becomes too sharp near the di- 
urnal maximum. The curves shown in Figure 15 of cos*' (4^/2) for various 
values of n show that the changes in the shape of the curve are large for 
small changes in n when 2 n ^ 8 and the shape changes very Jittle with 
n when n > 8 • The point to be made here is that when a powered cosine 
function is to be used to describe the diurnal density variation* the ex- 
ponent should be limited to values between 2 and 8. Indeed* Jacchia has 
consistently arrived at exponents in this range [5, 8* 10* 11, 12]. 

The U.S.S.R. Model is also a relatively simple, straightforward 
algorithm. Ir is very fast and relatively sophisticated; however* its use 

is limited to the altitude band from 140 km to 500 km and to solar flux 

-22 -2 * -1 

levels from 65 to 165 x 10 Watts m Hz . Furthermore, certain conditions 

can cause the model to yield negative densities. These conditions* which are 

physically realizable, occur when the 6-month average of the daily F^^^ ^ 

-22 -2 -1 

is near enough to 150 x 10 W ra Hz that F^ • 150 is chosen as the refer- 
ence flux. The scale factor which corrects the density for short term 

fluctuations in solar activity becomes negative for values of the daily 
and altitudes below the curve shown in Figure 16. It is true that the conditions 

for which becomes non-positive are not likely to occur often* but variations 

-22 -2 -1 

in - of the magnitude of 75 x 10 W m Hz have occurred and the 

potential user should be aware of this limitation In the model. The other 
factors are always positive in the altitude region between 140 and 500 km. 
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The Jacchia-Roberts models are the moat aophlaticated of the modala 
considered and provide upper atmosphere density va^ea over the greatest 
range of altitudea. The Jacchia-Roberta models and the U«S.S«R« Model pro* 
vide corrections for the diurnal variation » variations in solar activity 

O 

over both the ll*year and 27-day cycles » semiannual variation, and for 
variations In geomagnetic activity. In addition, the Jacchla-Roberts models 
account for seasonal-latitudinal variations in the assumed constant boundary 
condition at 90 km. and seasonal latitudinal variations in helium concentra- 
tions. The Modified and Asymmetric Modified Harris-Friester Model accounts 
for only the diurnal variation and the variation in density due to the vari- 
ation of solar activity over the 11-ycar solar cycle. It sould be remembered 
that the original model given by Harris and Priester [7] does include pro- 
cedures for accounting for all of the variations discussed herein. However, 
one of the fundamental motives in this investigation was the determination 
of an efficient model which can be adapted to model "real-time" atmospheric 
density variations and the Asymmetric Modified Harris-Friester Model has the 
desired characteristics of computational efficiency and adequate fidelity in 
representing the diurnal variation to form the base for such an adaptive 

model. The details of the formulation of such an adaptive model are given 
in Reference (6), 





7. Conclusions 


Three widely used atmospheric density models have been discussed 
The computationi^l aspects of each model have been shown and comparisons 
of the computitional speeds and computer storage requirements have 
been made* In general » all of these models can be said to be quasi- 
dynamic representations of the atmospheric density; that is» they are 
neither ruinpIcLo'ly sL.ii n* iiot: complotely dynamic* The time dopciulout 
variations in the model density profiles are determined by both the evalua- 
tion of explicit continuous functions of time and by the input of time vary- 
ing parameters to the algorithms. These input parameters » specifically 
measured values of solar and/or geomagnetic activity » arc available for use 
by the algorithms in discrete form only. Solar activity la reported as 
1-day averages and geomagnetic activity is available every 3 hours. Since 
the time delay between the measurement of a change in geomagnetic activity 
and the corresponding response in the atmospheric density approximately 
6.7 hours [5]» a direct data link with the geomagnetic activity index re- 
porting agency muld be required for real time or near real time applica- 
tions. Usually^ though^ some predicted average values are used with the 
ensuing errors being accepted as unavoidable* However » even if the input 
parameters are available within the lag time interval » the density models 
are still static with respect to the time interval between successive re- 
evaluations of the parameters and this lack of fidelity constitutes an 
error source in the evaluation of the drag forces* 

A method/ to overcome this shortcoming in the current density models 
is to provide a continuous input of measured solar and geomagnetic activity 
indices* However* such a solution would be difficult to implement In a neat real 



time mode. Another method to possibly accomplish accurate drag determine*- 
tions, especially in real time or near real time, is to estimate the drag 
by processing satellite observations with a sequential linear filter. 

This latter concept is attractive for a number of reasons. Besides allowing 
for real time determinations, the technique could allow for improved time 
an spatial resolution in the drag model, improve the performance of the 
filter by minimizing errors in the drag model, and significantly reduce the 
requirement for external data input. 
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Table 3» ^Sample Density Profile Determined bgr the 
Analytic Jacchia>Roberts Model Using 
Single Precision Arithmetic 


ALTITUDE 

(HETLRS) 

DENSITY 

m/v? X lo” 

300000.01 

2.16516 

300001.53 

2.16338 

300006.65 

2.16673 

300015.38 

2.16670 

300027.70 

2.16165 

300043.60 

2.15915 

300063.09 

2.15988 

300086.14 

2.16069 




Table 4. Comparison of Density Caloulattons with tbe 
Analytic Jacchia-Roberts Model 


ALTITUDE 

(KM) 

SINGLE 

PRECISION 

DENSITY DIGITS 
DOUBLE PRECISION 
EQUATION (A.12h) 
SINGLE PRECISION 
ALL OTHER CALCULATIONS 

DOUBLE 

PRECISION 

90 

3.46 

3.46 

3.46 

100 

5.4952 

5.4956 

5.4977 

150 

2.5798 

2.5800 

2.5809 

200 

3.9305 

3.9308 

3.9322 

400 

8.0247 

8.0253 

8.0283 

600 

4.2362 

4.2365 

4.2331 

800 

3.8025 

3.8028 

3.8042 

1000 

8.8327 

8.8333 

8.8366 

1500 

1.5922 

1.5923 

1.5929 

CENTRAL 
PROCESSOR 
TIME (SEC) 

.289 

.292 

.376 

CENTRAL 
MEMORY CORE 
STORAGE (WORDS) 

12000g 

12700g 

15300g 
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Table 6. Comparison of Densities Calculated by ttie Analytic 
Jaoohia-Roberts Model and tbe Modified Analytic 
Jacchia>Robeits for » 1100*^K 


ALTITUDE 

(KM) 

ATMOSPHERIC DENSITY DIGITS 

PERCENT 

ERROR 

UNMODIFIED HODFL 

MODIFIED MODEL 

100 

5.4977423 

5.4977547 

.0002 

no 

9.9303006 

9.9303229 

.0002 

120 

2.4596339 

2.4596394 

.0002 

125 

1.4018303 

1.4018202 

.0007 

200 

2.9381290 

2.9379562 

.0058 

300 

2.7866646 

2.7863754 

.0104 

400 

4.8761861 

4.8755852 

.0123 

500 

1.0416292 

1 .0414961 

.0128 

750 

3.6213252 

3.6210061 

.0088 

1000 

4.4213508 

4.4214982 

.0033 

1500 

7.6597326 

7.6603699 

.n083 
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COMPUTE 
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COMPUTE 

T(h") 

d.(h") 


DEFINITIONS 

T^ - exospheric temperatut'? 

Tj^ - atmospheric temperature at 12^^11 

h - satellite altitude (kml 

T(h) - exospheric temperature at h 

M(h) - mean molecular mass at h 

d^(h)- number density of i^" component 
divided by Avogadro's Number at h 

p(h) - mass density at h 


'h" < 125} 


COMPUTE 

T(h"') 

d.(h-) 

p(h'") 


Figure Lc^ic Flow Chart for Analytic Jacchia-Roberts Model 


63 








RADIUS (ICM) ALTITUDE (KM) DENSITY X 10 



6680 



TIME (SECONDS) 


60 ^ 100 


r 


1 

I 

r 


24 


« 


“r 

60 


-142 -113 -IT ~4 

LONGITUDINAL DISPLACEMENT (lEGREES) 

91 60. 30 0 30 45 30 

LATITUDINAL DISPLACEMENT (DEGREES) 
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Figure 6. Density Profile Variations 
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Figure 9* Density Profile Variations due to 

Changes in Geomagnetic Activity Index 
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Figure 12. Comparison of Modified Anal 3 rtic 
Jacchia-Roberts Densi^' Profile 
with Ai^mmetric Modified Harris- 
Priester Density Profile in a 
Nearly Circular Orbit 
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Figure 13, Comparison of Jacchia-Roberts (a) and 
Russian (a) Density Profiles at low 
Solar Activity Level 
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Atmospheric Density - KG/M^x 10 



Figure 14. Comparison of Jacchia-Roberts (a) 
and Russian (a) Density Profiles 
at Medium Solar Activity Level 
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A DEMONSTRATION OF THE VALUE OF 
GENERAL PURPOSE, ON-BOARD 
SATELLITE COMPUTERS 
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Space Analysis and Computations Group 
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I. INTRODUCTION 


The TRANSIT Improvement Program (TIP) satellites were designed 
and built by the Johns Hopkins University Applied Physics Laboratory for 
the U. S. Navy. These are navigation satellites which have onboard a 
general purpose mini-computer with 32 K words of memory. Also, each TIP 
satellite has a hydrazine- fueled Orbit Adjustment System (OATS), and an 
attitude control system which operates in both a gravity gradient and 
spin stabilized mode. The spacecraft is spin stabilized during the orbit 
adjust phase, and, later, operates in the gravity- grad lent mcde as a 
drag- free satellite. A similar drag compensation system (DliiCOS) was 
flcwn on the first satellite of the series, TRIAD. 

A picture of the fully deployed spacecraft Is shown In Fig. 1. 
During the initial orbit adjust phase, the scissors boom is folded, and 
the hydrazine rocket and tank are attached to the spacecraft. The 
four solar panels provide a configuration for stable spin about the 
longitudinal axis, lab led z in the figure. Later, after the hydrazine 
is used up, the boom is extended with the empty rocket system acting 
as an end mass for gravity- grad lent stabilization. 

The solar panels are designed to unfold Innediate ly after the 
spacecraft achieves orbit. When the panels failed to erect, the TIP-II 
spacecraft was left in a low power condition and with unfavorable moment- 
of-inertia ratios for spin stabilization. One year later, TIP-III 
experienced an identical failure. In addition, a boom deployment problem 
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later caused Che aclsaora boom links Co break on TIP-II under normal 
motor driven deployment. 

Under ordinary c Ire urns tancea , with herd-wire spacecraft logic, 
these problems would have precluded any parts of Che mission being 
achieved, and would have even prevented important engineering checkout 
of many of the on-board subsystems. However, the ability to change the 
flight computer software after launch allowed us to Implement various 
complicated work-arounds and achieve a partial mission success. 


4c 

♦ 



Figure 1, TIP-II Orbitai ConilguratlOQ 
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This paper describes how the fllf^t computer was quickly 
reprogrammed to perform various control functions which: 

1) performed power management to avoid troublesome space- 
craft blackouts; 

2) achieved enough spin stability to fire the OATS thruster; 

3) raised the parking orbit to a workable altitude; 

4) removed a high (45 rpm) tumble rate which was the Indirect 
result of one of the failures; and 

5) deployed the gravity- gradient boom successfully on TIP- III. 
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II. 


THE PLIGHT COMPOTER AND ITS SOFTWARE 


The TIP flight computer shown in Fig. 2 is a general purpose 
mini- computer with specialised Input/Output (I/O) logic to service 
various spacecraft functions in real-time. The computer consists of 
two redundant CPU's complete with I/O logic and two magnetic core mem- 
ories. Either memory or both may be used with either of the redundant 
CPU's. Each memory provides programmable storage of 16,384 words of 
16 bits each. There is also a 64-word hard wired, Read-Only Memory 
(ROM) containing a special loader program for restarting the software 
(Ref. 1-2). 


The TIP computer was desigxied for assembly language program- 
ming. The memory cycle time for the computer is 4.8 (jisec, with the 
time for an ADD operation being 9.6 usee. The TIP interrupt system Is 
a hard-wired priority system containing 32 inputs. The 24 highest 
priority intertupts are labeled external and the last eight are Internal 
As Implied by their name, external Interrupts are driven by systems 
external to and Independent of the computer. The eight internal inter- 
rupts are controlled by the software and are used for high-speed link- 
age to various subroutines. These interrupts can also be masked and 
enabled via software. 


All computer input data are transmitted via RF link. The 
satellite can receive digital data at a rate of 10 bps or 1000 bps. 

The slow rate can feed the computer ^ the command system, while the 
1000 bps data can only be used by the flight computer. There are a 
number of ways, direct and Indirect, in which computer outputs can be 
realized. Direct outputs occur when data from memory are transmitted 
by the RF downlink channel to the ground. Indirect outputs are inferred 
when another satellite sub-system changes in response to a directive from 
the flight computer. The most useful direct output occurs In the com- 
puter dump mode. Upon command, the telemetry system transmits contin- 
uous flight computer data (via TM modulation) . This mode requires a 
dump program In the computer to relay the contents of mei>iory to th* 

TM system at the proper rate. 


82 






ORIGINAL PAGE IS 
DE POOR QUALITY 



Figure 2. TIP-n, Computer Proceasor, S/N TP-1 
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Th« fliiJtt computer loftware conelets of « lyetem of inter- 
rupt driven* real-time prograa». These programa perform on-board data 
nagkment and Interact through special hardware Interfaces with other 
subsystems to give the computer far-reaching powers in controlling and 
monitoring the aatellite. 

The flight computer derives most of its power to perform 
control functions by virtue of its direct interface to the spacecraft 
telemetry and command systems. Am rXF Mfdware includes a telemetry 
(TM) system whose function is to gather* precMS and format space- 
craft data for transmission to the ground in a serial bit stream. 

The TM system Is digital* with $ bits per chanael* 172 channels per 
frame* and a 4.227 sec frane rate. The TM Interface allows the com- 
puter to exchange data with the TM system under direct software control. 

To receive TM data* the software requests via the interface 
one of the 172 TM addresses. When this address <^curs within the normal 
cycling of a TM frame (every 4.3 sec)* the computer is interrupted and 
receives the data for storage or processing. Generally* data stored in 
the flight computer mei<*-ry is later returned to the in the 

form of a memory dump transmission. 

The TIP command subsystem contains digital (10 bps) logic to 
perform ;:he remote execution of relay commands, pulse coomiands* digital 
data commands and slow (10 bps) loading of the computer memories. Through 
the comnand interface* the flight computer has direct access to the front 
end of the command system. Any connand can be issued by the flight soft- 
wara by serially transmitting the comnand bits through the interface at 
the required 10 bps rate. The length of a relay command bit string 
requires 2.3 secs for complete transmission. Any comnand can be executed 
with a progx’amned time delay by allowing the computer to issue the com- 
mand. This "delayed comnand" capability results from loading the informa- 
tion for the delayed commands into the computer memory to be processed 
at a programmed time. 
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The main implications of the I/O Interface described above 
are fi. ; the computer la limited to a data sampling rate of 4.3 s»c 
for any given TM channel, and the maximum command rate la one every 
2.4 sec. These conrtralnta became quite important 'in some of the con- 
trol functions Implemented. 

Tlie TIP ground support system Is illustrated in Fig. 3. The 
software for this system Includes at least five major programs and 
uses four different computers. An overview of the ground system Is 
contained In Ref. 3. The back1>one of the system Is the ground station 
PDP-il-40, operating through a front end PDP>8. This system is used to 
control all real-time satell'tte operations, and Is also used for data 
formatting, real-tliue conversions and display, and miscellaneous ucilltins. 

A program to be Injected Into the TIP flight computer begins 
as a card deck which contains the program code written In the flight 
computer assembly language. The card deck 1^ Input to the IBM 360/91 
computer and processed by an assembler program called ARTIC. The output 
of ARTIC Is the machine cede on a magnetic tape along with a printed 
listing of both the Input assembly language Instructions and the corres- 
ponding machine code. The program tape Is then stored on a disk file 
In the PDP-11 by the TIPLIB program. This disk file library contains 
the latest versions of all the flight computer software, Including 
operational and diagnostic programs. 

The FdP-ll program that selects flight computer programs from 
the library and formats them for transmission to the satellite Is called 
TIPLOAD. The Input to TIPLOAT Is a card deck which defines the pro- 
grams to be selected from the disk file library. This data from the 
library is then mergid with other flight operation Inputs and formatted 
for transmission to the spacecraft. The output of TIPLOAD Is a disc 
fll.5 (LDM file) In the PDP-11. The data on this file Is arranged Into 
segments called "modules" which can later be Individually transmitted 
to the spacecraft. 
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Operations 



Figure 3. TIP Ground Support System 














During a satellite pass the L!M file data Is transmitted to 
the TIP spacecraft under control of the TIPCOM program, which also 
resides In the PDP-ll. In addition to transmitting data up to the 
flight computer, TIPCON also receives and records downlink loading 
feedback flags from the flight computer. All real~time communications 
are handled through TIPC(W. In addition, TIPCCM converts and displays 
on CRT much of the normal TM data In real time. 

The overall ground system Is complicated, but very flexible. 

It gave us the ability to completely reprogram the flight software after 
launch, as well as to manage the system In orbit In ways we had never 
dreamed of %dien the software was developed. 

The main flight computer software Is a set of basic programs 
called SYS which are resident In memory at all times. SYS contains 

1. loading programs which can handle data at 10 bps 
or 1000 bps; 

2 . a memory dump progr A which can read out areas of 
memory on either a 325 bps or 1300 bps downlink; 

3. a status routine which sends 80 bits of computer 
Information to the TM system each TM frame; 

4. a timekeeping routine that keeps a high precision 
universal time (UT) clock. The basic unit of time 
Is referred to as a "tock" and Is precisely 120/6103 
seconds ; 

5. a Thne-Queue program which controls the chronological 
sequencing of computer events such as delayed commands. 

For details on the complete flight software system, see Ref. 4-5 

In addition to SYS, there are other special programs which are 
loaded by SYS when needed. TWo of those programs are 

1. Delayed Coimaand Program (DCPRO) 

'*Delayed cosaMiid" refers to a relay or data casamnd which Is 
sent by the flight computer directly to the satellite command system at 
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socH prc~«peclfled tloe. The delayed eooBaiade are prepared by special 
card Inputs to the TIFIOAD Processing Prograsi. The TIPLOAD Processor 
autonatlcally fonaits the coonand bit strings and the appropriate Tine 
Queue entries as needed by the fll^t computer. Note that the coonand 
system hardware requires that there be at least .1 seconds between com> 
mends. DCPRO may also be used to send commands upon the occurrence of 
certain eventa. This Is done In conjunction with the next program to be 
described, TMDN. TMON Initiates delayed coonands whenever the data In 
certain telemetry words matches prespecified criteria. Since It takes 
2.3 seconds to send a relay command, the fll^c computer ttoftware must 
make sure that the extat-nsllv triggered delayed commands do not Inter- 
fere each other or with the time ordered ones being controlled by 
the Time Queue. 

2* Telemetry Storage Program (TMON) 

TMON Is used to sample and store real-time telemetry data In 
Che flight computer. The program expects as Inputs: 

1) A start time (Time Queue Entry); 

11) A list of TM channels to be stored and the rate at 
which each Is to be stored. 

TMON allows each TM channel to be sampled at Its own rate, hence all 
channels need not be sampled during the same frame. 

The program remains operating In the flight computer until the 
specified storage area fills up. The program automatically stops storing 
data at this point. Once the program has completely executed, a memory 
dump procedure Is needed to transfer the stored data from the flight 
computer to the ground station. 

The use of the Time-Queue for delayed commands, and the ability 
to send coiunsnds while monitoring TH functions proved to be extremely 
valuable after the TIP failures. In our wildest imagination we could not 
have foreseen the use we would make of these programs, nor the salvation 
they would provide for the crippled mission. 
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III. 


THE EFFECT OF THE FAILURES 


The TIP spacecraft are lar.nched by Scout Vehicles Into a 
nominal polar, 180 x 400 n. ml. altitude parking orbit. The orbit adjust 
system Is then used to change tbA orbit to circular at 600 nm altitude. 

At the same time the inclination Is precisely trimmed to a selected 
value near 90° to control the nodal precession. An important part of 
this operation Is to select optimum directions for thrusting to cor> 
rect the altitude and Inclination together and minimise the fuel require- 
ments. 

The spacecraft Is dealgmd to be spin stabilized about Its 
longitudinal symmetry axis (z axis) to provide stable directional con- 
trol during firing end to compensate for thruster misalignments. To 
achieve spin, an analogue magnetic dipole spln-up system provides con- 
tinuous torque about the z axis using the earth's magnetic field. 

Passive nutation dampers on the end of the solar panels negate the 
effect of random transverse torques introduced by the spln-up system 
during this operation. To slew the z axis to a desired firing direction, 
once spin Is achieved, a reversible z-dlpole coll Is aboard to provide 
precesslonal torques using the earth's magnetic field. 

After the orbit Is adjusted, the remaining hydrazine Is 
vented, and the empty OATS system becomes the end mass on a scissors type 
boom for gravity- gradient stabilization (see Fig. 1), 

When the solar panels failed to deploy on TIP-II we were left 
In the following situation: 

a) the spacecraft was generating less than half normal 
power; 

b) the spacecraft was not stable In spin about z; 

c) the nutation dampers were not In the correct position 
to be effective In dsmplng out torques transverse to 
a; and 
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d) the 60 lbs of liquid hydrazine at the top of the space- 
craft was an effective nechanisn to quickly transfer 
any spin about z into tusd>le about the stable trans- 
verse axis. 

For various reasons, the spacecraft needed to he at an average 
altitude of at least 400 n. tai. to be able to operate effectively as a navi- 
gation satellite. Also it was necessary to reduce the eccentricity to 
achieve the gravity-gradient stability required to do a good DISCOS 
experiaient. 

As will be described, ve managed to achieve this, leaving 
about half the hydrazine to be vented before putting out the boom. 
Unfortunately, the hydrazine venting system was designed for stable 
z-axis spin; and with our configuration a tumble torque was inevitable. 

We had no way of knowing how bad this would be, and were forced to 
take our chances and vent. The next time we saw the spacecraft, the 
tumble* rate was 45 rpm and the solar panels had been ripped free by the 
centrifugual force. This improved the power situation, but before the 
boom could be deployed, the tumble motion had to be dissipated. This 
was done by implementing an interesting digital phase- locked- loop and 
using the spacecraft z-coll to work against the earth's magnetic field. 

When we finally attempted to run out the boom on TIP-II, the 
links broke because of an unforeseen problem with scissors booms. 
Unfortunately, this happened after TIP-III was launched so the problem 
had not been corrected. It was possible, however, to work around the 
problem on TIP-III by using centrifugal force generated by a tumble 
motion. W^ were able to solve this problem by reversing the de- tumble 
program mentioned above. 


* 


Throu^out, we will refer to tootion about the spacecraft longitudinal 
(z) axis as "spin," and motion about a transverse axis as "tuad>le." 



IV 


POST UUNCH OPERATIONS 


1. Power Management 

The Imnedlate problem after the solar panels failed to deploy 
on TIP-II was the severely restricted power capability of the spacecraft. 
This problem was particularly bad when we tried to use the magnet:*c 
system for spin-up and precession, and the problem was exacerbated by 
the fact that we launched into a minimum sun orbit. 

To protect the battery, the power system is equipped with a 
low-voltage sensing switch, (LVSS) which shuts down the main power bus 
automatically when the battery voltage reaches 13.8 volts. Although 
this is not a disaster, it was extremely inconvenient when it occurred 
because the spacecraft had to be restored to its previous state through 
a complicated series of commands. The spacecraft ^dilator shifted 
frequency dramatically during the first few minutes of warm-up making 
it very hard to keep the receivers locked up. Also, the flight computer 
system had to be re- loaded, restated, and the U. T. clock reset: each 
time the LVSS tripped. With passes only about 8 minutes long, this 
led to a hectic operation with the spacecraft frequently rising silent 
and not responding immediately to the recovery commands. 

To solve the problem, we modified the TM storage pxogram to 
monitor the battery voltage channel once per minute. When the voltage 
fell to a threshold level, the program used delayed commands to throw 
off the magnetic system power. The threshold was inputtable ani# usually 
set to 14.5 volts. The magnetic system draws about 60 - 70 watts, and 
relieving this load when the battery got to 14.5 volts was generally 
sufficient to prevent the LVSS from tripping. 

Even with the battery voltage monitored, the spacecraft systems 
had te be duty cycled to prevent power drain. The {ime-Queue feature 
of the flight software, discussed previously, was used to iturn systems 
off and on at scheduled times to effect the duty-eyeing. Hogevai:^ 



this proved to be a great deal of work* -punching cards « preparing 
updated Time-Queue files with the TIPLOAD program, and Injecting the 
data into the spacecraft memory. When the Time- Queue software was 
developed. Its primary purpose was to fire the OATS rocket out of view 
of a ground atatlon, and It was not designed to easily handle a large 
number of delayed comnands (like 100/day). The operations team quickly 
found most of their time taken up trying to keep the flight computer 
fed with duty cycle data. We were really spinning our wheels. 

It proved rather easy to reduce this workload dramatically. 

We quickly recognised that nearly all of the chrono! ^glcal delayed 
commands for duty-cycling were periodic in time. They were generally 
tied to the orbit geometry, such ns: systems turned off and on over 

the equator or the poles, systems turned off In the earth* s shadow, etc. 

The answer was to make the Time-Queue automatically cycle Itself. 

A program called CYCI£ was written to cause the Time-Queue 
actions to repeat with an Inputtable period. This Is done by calling 
CYC1£ via the last Time-Queue entry in the list. The CYCLE program sets 
approprla«.e pointers to the stsrtlng conditions and then restores the 
original Time-Queue list, adding the input period to the time for each 
entry. This causes the Time-Queue actions to be periodically repeated 
Indefinitely, until the process Is stopped by ground command. 

• 

Some people \rare nervous at first about relinquishing control 
to the computer and allowing the spacecraft to operate autonomously. 
However, this simple fix worked beautifully and they soon became believers. 
Even the simplest operation would have been difficult to carry out without 
CYCLE; and, as will be seen. It would have been next to impossible to 
carry out the more complicated operations we eventually undertook. 



2 . 


Spin-Up Operations 

The early operations with TIP-II were Involved with attempts 
to spin the spacecraft about the z-axls. There were two reasons for 
these attempts: (1) early In the game we had hopes that sufficient 

spin could generate enough centrifugal force to break the solar panels 
loose; and (2) stable spin was required to fire the OATS rocket. 

The spln-up system Is a feedback control system that uses the 
earth's magnetic field to torque the spacecraft. Two orthogonal colls 
(x and y) provide a dipole moment In the spacecraft to supply the torque. 
The earth's field Is continuously sensed, and the x and y coll currents 
are automatically phased to keep the resultant dipole orthogonal to che 
component of the earth's field that lies In the x-y plane. (The geometry 
Is shown In Fig. 4.) This results In a torque about z that Is always In 
the same sense, along with random torques transverse to the z axis. 

With a stable configuration the spin about z will gradually build up to 
the desired level, while pass^e dampers remove the nutation Induced 
by the transverse torques. 

With our unstable configuration, It was a different story. 

When the spin system was turned on, the liquid hydrazine was sloshed by 
the transverse torques and acted as an effective mechanism to quickly 
transfer any spin into tumble motion. Initial attempts to acl^eve any 
spin above a 1/4 rpm were unsuccessful. To help the problem, we modi- 
fied the^fllght software to control the times that the spin ^stem was 
^n, so as to minimize the transverse torques. 

The spacecraft are equipped with 3 orthogonal magnetometers 
measuring the body-fixed x, y, and z components of the earth's field. 

We modified the TM storage program to sample the channel Upr the z com- 
ponent every frame, and use It as criterion for tilrnlng on the spin 
system. Each TM frame the program made the following test: 

• |Mj < c 

Where Is the z magnetometer read^g, and 

c Is an Inputtable threshold. 


Spacecraft Longitudinal Axle 
z 



Figure 4. Magnetic Spln-Up Geometry 


When the threshold was satisfied, the program used a delayed cooenand to 
turn on the spin system, and conversely turned the system off when the 
z component was out of range. We had to Include additional logic to 
prevent the program from continuously sending comnands once the system 
was In the correct state. By allowing the system to be on only when the 
z component was near zero, the transverse torques were nearly eliminated. 

It happens that there are two relays In series to control the 
magnetic system power. Ihls turned out to be quite useful since we 
could let one relay (normally on) be controlled by the battery voltage 
and the other be controlled by the z magnetometer reading. Thus we 



could allow both monitoring functions to operate In parallel, keeping 
our protection against the LVSS tripping during spln-up. This turned 
out to be about the only piece of good luck that we had In the opera- 
tion. 

It took about two days to generate the spln-up program, and 
after some trial and error we settled on a value for c of about 10 
per cent of the full-scale field reading. The program helped quite a 
bit. We were able to achieve spin rates up to 1 rpm, but then the 
spacecraft would gradually build up nutations and transfer to tumble 
with a time constant of about 1 orbit (90 minutes). This was still not 
enough spin rate or stability to fire the OATS rocket. 

At this point the Idea arose of using the spacecraft z-coll 
as a device for actively damping the nutatlonal motion to maintain more 
stable spin. This Is the coll that Is normally usad to provide pre- 
cesslonal torques for pointing the rocket nozzle. The coll can be 
switched by relay command to either the plus (dipole along +z) or the 
minus (dipole along -z) state. 

For clean spin without nutation, the spacecraft magtietometers 
record a very distinctive pattern from the earth's field vector. Since 
the earth field Is varying rather slowly due to the orbital motion, the 
attitude dynamics dominates the magnetometer's variations. The x and y 
colls record a sinusoidal variation 90° out of phase with period equal 
to the spin period, and the z magnetometer records a nearly constant 
. value. As nutation (coning) builds up, the z reading begins to show 
an oscillation with period equal to the nutation or coning period. 

The Idea behind the damping program (DAMP) was to let the 
computer sense the derivative of the z magnet'sineter reading, and then 
set the polarity of the z-coll to produce a transverse torque (using the 
earth's field) that opposed the derivative. Again we '‘could modify the TM 
storage program, this time to control the z-coll polarity based on the z 
magnetometer reading. The logic was simple. 



1. Bach frame, compute the difference between the 
current z magnetometer reading and the previous 
frame's reading, as an estimate of the derivative. 

2. When the difference changed sign from plus to minus, 
send a delayed conmand to set the z>coil polarity to 
a minus dipole. 

3. When the difference chatiged sign from minus to plus, 
send the command to change the z>coil to a plus dipole. 

We also had to include additional logic to correct the 131 reading for 
the z«coil effects. The z-coil strength is of the same order as the 
earth's field at 800 km altitude so It has a large effect on the z 
magnetometer reading. 

Note also that the DAMP program can be used to de> tumble the 
spacecraft when it is in stable tumble about a transverse axis. In this 
case the z magnetometer records a sinusoidal variation. The effect of the 
program is to continuously adjust the z>coll polarity to produce a torque 
opposite the motion. This proved to be quite handy in reducing the time 
required to dissipate tumble motion. 

The magnetic system is designed so that the spin system or the 
z*coil can be in use, but not both modes at the same time. As soon as we 
had the DAMP program written, we quickly found that it worked well, but 
we needed to continuously alternate between the spin-up mode and the 
damping mode to be effective. This meant we needed a way of dynamically 
switching between the SPINUP logic and the DAMP logic in the TM monitor- 
ing program. 

The idea of changing the program logic oynamically while the 
flight software was actively running was a completely new feature for 
our software system. It was not hard to implement, and it turns out to 
be a very powerful capability. 

We compiled both sets of logic (SPINUP and DAMP) into the TM 
monitoring program, with a simple program switch to select which path 
would be used to evaluate the z magnetometer reading. We then wrote a 
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program (ALTFRO) i^lch was called by a Tine-Queue entry to change the 
switch at a scheduled time. ALTPRO was written to be rather general, 
driven by an Input list of nenory addresses and their corresponding new 
contents. Each time the program Is called at a scheduled' time. It 
works on the next entry In the list, putting the new contents Into the 
specified addresses. 

To accomplish a splnup, then, the scenario ended up something 
like the following. At a pre-scheduled tine, with the spacecraft 
essentially motionless, the computer would select the magnetic system 
for splnup, turn on the TM system, and activate the SPINUP logic. 

Near each equator, the computer would select the z-coll by delayed com- 
mand and switch to the DAMP logic using the ALTPRO program. After 
several minutes of damping, the spin system would be re- selected and 
the logic switched back to the SPINUP program. This entire set of 
Time- Queue entries were then cycled with the orbital period by the CYC IP 
program for continuous operation. Of course, during the complete scenario, 
the battery voltage was monitored to prevent the LVSS tripping. If the 
battery monitoring program did shut down the magnetic system, the 
Time-Queue cycle would turn It back on at the beginning of the next 
cycle. 

By starting with a fully charged battery, the above scenario 
achieved spin rates up to 4 rpm, which was enough to be able to fire 
the OATS rocket. However, we still did not have directional control 
over the spin axis. If we left the spacecraft alone after achieving 
3-4 rpm, it would maintain Its spin stability reasonable well for about 
one orbit (90 minutes). After that It would begin coning, and once It 
started. It would go quickly Into tumble. As soon as we tried to process 
the spin axis with the s coll, the Induced nutations would make the 
motion become unstable. 

We tried alternately processing and damping, but we rapidly 
ran out of power*. There was no way to get directional control, and 


and also spin. The damping program continuously robs the spin of 
energy. 
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without It the orbit edjustnent looked very bleak. The solution to this 
problen required ingenuity, hard work, and nerves of steel. 

3. Firing the OATS Rocket 

The basis for the tine and inertial direction of the OATS 
firings is n solution to an optimization problem to minimize fuel for 
the required orbit change. Our normal procedure 's to solve for the 
true anomaly and the direction of thrust from the current orbit and desired 
orbit parameters. This gives us three degrees of freedom to solve for 
each firing, and thereby optimally correct the semi-major axis, the 
eccentricity, and the inclination. Without directional control we could 
not hope to use this scheme. 

After some thought, however, we came up with a "far-out" idea 
that we suspected might work. Each time we ran the spin- up scei lo we 
cane up with a distinctly different inertial attitude, and although we 
could not change it, the spin axis direction remained gyroscopically 
stabilized for about one orbit. We could make use of this fact in the 
following way. 

Instead of solving for a time and direction of thrust from the 
orbit parameters, we could accept the direction we had, and solve for 
the optimum time to fire, given that attitude and the orbit. We could 
then obtain a measure of how effective the thrust would be by comparing 
the resulting orbit changes to those we would have achieved if we could 
have chosen the direction. If this measure was reasonably high, we 
would fire the rocket. Otherwise we would de- tumble the spacecraft and 
try again. With the spacecraft coming up in a random attitude each time 
we spun it up, a certain percentage of the time we would get lucky and 
obtain a favorable attitude. 

Several problems needed to be solved before the above scheme 
could be implemented. First of all, the calculations for the optimum 
firing point had to be done in real time during a single pass because 
the spacecraft would not spin stably longer than about one orbit. 
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Thus, He deeded to begin an automatic spinup scenario four or five hours 
before our pass, such that the spacecraft rose spin-stabilized at 3 - 4 
rpm. Then In real time during the pass «e would: 

1. Determine the inertial attitude of the spin axis. 

/ 

2. Calculate from this and the current orbit the optlzHim tine to fire 
In the next 90 minutes. 

3. Make the decision whether the thrust would be effective enough. 

4. Inject the appropriate Tine- Queue and Delayed Conniand data Into the 
spacecraft to control the firing at tdie specified time. 

We had about 8-10 minutes during the pass to accomplish the above 
operation. 

Luckily the first requirement was already satisfied by an 
existing capability. The TIP attitude determination software had been 
designed so that It could be run In real time during the satellite 
passes. We begin by describing this system. 

The attitude calculation (and the orbit calculation as well) 

Is too complex to be handled by the PDP-11 ground station computer. 

This computer has Its hands full during the pass handling the satellite 
data link. The attitude calculation Is done on the large IBMh370 com- 
puter which Is connected by telephone data link to the PDP-11. The TM 
data for the attitude Is fed In real time to the IBM-370, where the cal- 
culation Is done Interactively In a "time-shared" operation (TSO) se8Si.on. 
The system Is shown schematically In Figure 5, and you will note that 
we had the capability to operate throu^ a station In Hawaii as well as 
the AFL ground station. In this set up, the PDP-11 acts as a TSO terminal 
for the 370 computer, controllirtg a second TSO session which receives the 
raw TM data atid passes It throu^ a shared disc file to the attltude/orblt 
TSO session. 

To use the system for our scheme, we had to add to the 370 
attitude software the extra program to handle the orbit calculation des- 
cribed as item 2 above. The required equations are developed and dis- 
cussed In Appendix A. This turned out to be a non-trlvlal program. 
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The software was run Interactively and the Infomatlon about the firing 
was displayed on a CRT graphics temlnal In real tine. A sample of the 
output display Is shown in Figure 6. 



Figure 5. Real-Time Computing System Used for OATS Firings 
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OATS FtRXNQ IHFORdATXON 


THE niRN FXCURE OF HERXT • 
S-T-tt DIRECTXON AT FXRE • 


CHANGE TO SENIHIAJOR AKXS • 
CHANGE TO ECCENTRXCXTY 
CHANGE TO XNCLXNATXON 
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UT TXNE FXRXNQ 
TINE OF FXRXNG 
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50940 (SEC) 

76.00 cnXNUTES AFTER SET) 
836A«5F98 


HIT CARRXAGE RETURN TO RETURN 
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Figure 6, Real-Time CRT Display for OATS Firing 


From the program display, the decision was made whether t 
fire or not. We generally had 3-4 minutes to run the program and 
make this decision. The decision was based on many factors, not the 

least of which was whether we believed the answers we were getting from 
the software. 

At times the decision was difficult. For example, in the case 
where the optimum firing time was nearly one full orbit after set, this 
meant we had Just passed the optimum point. The question then became: 
should we wait for nearly 1-1/2 hours to get to the optimum point and 
risk the inevitable nutation build-up? Or, should we fire immediately 
during the pass while we had good stability and accept somewhat less 
than optimum geometry? This type decision had to be based on a multi- 
tude of factors such as the current spin rate and nutation angle, how 
far past the optimum point we were, and how good the thrust would be if 
we waited. 
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There ves one final problem to be solved before the entire 
scheme could be made to work. This Involved Item 4 above ~ transmitting 
the firing data back to the satellite computer In real time. As men- 
tioned in Section 2, the Time-Queue and Delayed Command data are prepared 
and formatted for transmission by the PDP-ll<jpr 5 ^am called TIPIXUU). 

The problem was that TIPLOAD was not designed to be run in real time 
during a pass. It was a pre-pass utility. And even If it could be run In 
real-time, the PDP-11 was completely taken up during the pass by the 
TIPCC»I program handling the real-time data links. The UPLOAD function 
Is a formidable task and there seemed to be no way to get In the Time- 
Queue data to control the firing. 

The answer turned out to be beautifully simple. We made all 
of the times In the calculation and In the Time-Queue be relative to 
the satellite set ^Ime for the pass. Then we simply controlled the 
firing time by setting the satellite clock to a dummy value to make the 
rocket fire at the correct amount of time after set. The clock was already 
designed to be easily set in real time by simple keyboard type-in. The 
Time-Queue and Delayed Comnand data for firing could be formatted as a 
fixed set of times, and prepared for Injection once and for all In a 
single TIPLOAD run. 

The actual OATS firing passes were quite complicated. A suc- 
cessful operation required that all of the computers Involved as well as 
the telephone data links be operable for the pass. We had Tnany ope^tlons 
scrubbed because of computers going down, or because of simple human error 
caused by the time pressures. A synopsis of the entire operation Is as 
follows : 

1. On the day prior to the pass, transmit a spln-up scenario Into 
the flight computer, scheduled to begin five hours before the 
pass. 

2. Just before the pass, bring up all computers, establish the tele- 
phone data links In full duplex, and Initialize the T5^0 sessions 
on the IBM- 370. 
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3. As soon as the satellite rose, send conmands to begin transmitting 
TM attitude data. At this point we had less than 10 minutes to com- 
plete the operation. 

4. Record attitude data on the IBM>370 for about three minutes and then 
begin the attitude and orbit computations. 

5. While the orbit computation is being carried out. Inject the 
"relatlvc-tlme" delayed command and Tine-Queue data Into the 
flight computer. 

6. Make the decisions about the firing. 

7. Either fire immediately, set the satellite clock to the correct 
dummy value, or abort. 

8. On the next pass, set up the computer to begin a two-day de- tumble 
operation using the DAMP program. 

The scheme worked well, although It was a trying experience. 

(On one harrowing pass we actually fired the ro-^ket backwards, but all 
other firings were successful.) We were able to average about two four- 
mlnui'e firings a week for a month or so, using up nearly half the fuel. 

But. then, tUi.ngs began to go badly. We had always felt that a half-empty 
fuel tank would cause worse stability problems than a full one simply 
because there would be more sloshing around of the hydrazine. Sure enough, 
this oegan to happen, and worse yet the spacecraft began to come up spinning, 
consistently oriented normal to the orbit plane. We attributed this to 
spln-orblt comilng which got worse with the Increased damping of a half 
empty tank. Since we were most Interested In raising the altitude, this 
geometry was totally unfavorable. After trying without success for several 
vee’'s, we realized tha^ a new idea was needed. The Idea was not long In 
h..^ming, and once again the on-board computer saved the day. 
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4. The Tumble-Thrust Program 

The new idea for firing the OATS rocket was a radical departure 
from our previous method - we completely abandoned the idea of trying to 
maintain any attitude control during the firing. Instead we decided to 
let the spacecraft ^pible, since tha!| is what it wanted to do, and fire 
the rocket in short bursts when it happened to be pointing in the best 
direction. 


Tims, we could let the flight computer continuously determine 
the spacecraft attitude and then quickly fire the rocket by delayed com- 
mand at the correct times. Attitude determination programs are non- trivial, 
and this would have been an extremely difficult program to write had it 
not been for some simplifying circumstances in our case. 


At that point in time we were getting close to the orbit we 
wanted; we needed only to raise the perigee altitude some more. Thus we 
were willing to do all our firings In the along-track direction In the 
vicinity of apogee (within 30'* or so in true anomaly). This would raise 
the orbit without increasing the eccentricity and without affecting the 
Inclination. It turns out for a near polar orbit that there is a simple 
relationship between the magnetic field in the equatorial regions and the 
along-track direction. This geometry is Illustrated in Fig. 7. 

It can be seen from Fig. 7 that the field lines are roughly 
parallel to the flight path in the regions near the equator, so that ^ 

when the spacecraft is aligned with the field it is aligned with the 
velocity vector. Hence the approximate determination of alo.ng- track 
orientation in these regions becomes trivial using the sampled magneto- 
meters. 


We had observed during our prior operations that, when tumbling, 
the spacecraft angular momentum vector tended to align Itself normal to 
the orbit plane. Apparently there was a strong spin-oibit coug|.ing caused 
by the perturbing torques. In this case the tumble motloig is in the orbit 
plane, and it continuously carries the longitudinal axis through the space- 
craft velocity vector giving ample along-track firing opportunities. 



Earth spin 
axis 


Magnetic field 
lines 



Equator 


Orbit path 


Figure 7, TIP Orbit Geometry for Tumble-Thrust Program 


To control the firings, we developed a program called THRUST 
which was another special version of the TM monitoring program. Ihls 
program continuously monitored all three magnetometer channels, and 
determined when the spacecraft was aligned with along- track by the follow- 
ing test: 

• 

|M 1 < C, and Im 1 < C. and KM < 0, 

' x' 1 ' y ' 2 z ' 

• 

where M , M , and M are the orthogonal body- fixed magnetic field 
X y z 

readings, and are Inputtable thresholds, and K ■ *fl for a 

north-going geometry or -1 for a south-going geometry. The test on M 

z 

Is to establish that the thrust direction will be parallel to the velocity 
vector, rather than anti-parallel. When the tests were satisfied the pro- 
gram would Imnedlately Issue delayed commands to fire the rocket for two 



• • 


second* and then rearm for another firing. The CYCLE program was •s^d 
to activate the THRIJSt logic periodically for about 10 minutes nea» each 
equator crossing closest to apogee. « 

The ^ogram worked best at a low tumble rate because af the 
TM sampling ri^e and the 2.3 second time required to issue the faring 

command. For example. If the tut^le rate were 1 rpm, the thrust axis 

. • • o d 

would move about 14 In the time needed To Issue the command ahd another 

12*^ during hhe i second burn. This would amount to an effective 20° 

0 

pointing error for the thrust. Also, at 1 rpm the spacecraft would move 

• about 26° between TM samples a)^ the 4.3 second frame rate. Thus we had 
to ^pen the thresholds C^ and C^ to a reasonably high value (about'^0^ 
of full scale) to Insure an angular window large enough to avoid Msslng 
opportunities due to the TM sampling rate. • 

• We seeded only approximate along-erack orientation (+ 45°) on 

each firing to do the job, relying on the off-track components to cancel 

* from firing to firing. The TM sampling times tended to fall randomly in 
the angular window established by C^ and 'C^ , and this helped average 
the off-track components. However the command time lag caused hisses. 
These were not too serious since the off-tr^k components tended to be 
radial thrusts which only affected the eccentricity a bit. 

^ Because of the above considerations, about one rpm was the * 
pracCItcal limit In tumble rate for the program to work effectively. We 
started the spacecraft tumbling very slowly, and found that the firings 
themselves caused the tumble rate to continuously Increase, due to the 
^ small displacement of the thruster axis with the spacecraft cenCeV of 
mass. The rate Increased about 1 rpm for each minute of fifing. Thtfl 
meant that the scenario had to periodically switch to the DAMP prograni^ 
to de* tumble a*d prevent the tumble rate from building up. 

f'’ 

,0nce %he scenario was worked out, the process worked well. 

We simply sat back while the flight ctSputer continuously pushed the 
altitude up for a week or so (The people determining the orbit were 
qMfte startled when the period began creeping up orbit by orbit. Their 




software could not handle that case.? ff em «ii taitial parklij orbit 
o^ about 180 perigee altitude bp 38Q ti.<mi. apogee altitude, <wl * 

finally achieved an orbit of about 320 n.ml. 450 n.mi. ^ • 

At that point, it was important ta vent the remaining* hydro* 
zlne so that the spacecraft could be put 4mto the gravity gradient mode. 

The next spacecraft, TIP>IIX, was scheduled i&t launch^ «nd we wanted to 

• • 

continue the engineering checkout so that otheer <poCcnC|at ^roblema <v«wid 

be uncovered before* it was launched. The weating eastern vas 4eatgned • 

• • • * 
te release the hydrer<ne to the side with the spacecraft spinning a^out 

the* 2 axis. Since the spacecraft was unstable in spin, a large tumble 

motion resulted from«the venting operation. . * 

• • . * • 

• * . The tumble .motion was much faster than we had. hoped for.. The * * 

•• 

solar pane la were wrenched free by the centrifugal force end the space- * 

craft ended up tumbling at 45 rpm. This was being dfsnlpated at the . 

• ♦ * • • 

rete of about .2 rpm per day because «f the spacecraft magnetic hysteresis. 

The gravity gradient boom could *not be cVected until the spaecrevaft was* 

• • 

stationary, and we needed to do this IH teae than three months to^pro* •* 

perly lead the TIP- III. launch. We were ii^an obvious time bind,, atpl .* ** * * 

^ * • » * 

needed to> use the z coll to help dissipate the tumble motion. The pto- 

* ♦ . 

‘ gram developed for this turned out to be the most complicated of git the 
. ***** 
special post* launch. programs. « * 

* 

* 

‘ • 
S. *. ’ A Digital .Phase- Locked Ikx>p for Oe- Tumble 

« • * 

Ihe DAMP program, which switches the z-coil to .produce a de- 
tumble torque,, begins to lose effectiveness over about 1.5 rpm. At 
2‘rpm it does not work very well, and at 45 rpm it is useless. Some 
. simple arithmetic shows the problem, tt takes *t^o TM frames to sense the 
derivative** in the z magnetotdeter,'and then another 2.3 seconds to 
send the command .to switch polarity. Ihe resulting average time delay 


• • 
$ 


• ^ 




« * 


# 

P 
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is 6.6 seconds^ and at 2 rpm %kls £Must|s fMf9* th# ^hasa Qf the 

^ ® A A# 

tumble motion. At 45 rpm^ flid period ed the miotlisa As f »m «e|bafls so 

that even continuous consandln/ cannot swlteh the a-col^ ^elavAhd ^ce 

• e 

per cycle. • e 

• • 

A program was needed* thaC could determine the phase «f the 
^b|e motiOi^ well*ftfieaghi VtiMtatlClpate the peaks and t^aughs In the ^ 

z magnetometer tepdlng. Command strings could be started with 

• • 

. ^ precisely %he lead time needed to have the conmands 4ake at the right 

• • • •• • 

• AnStanC. For motions faster ghdig Sf Apmuhere Ahe period Is less than 

the j|.4 second ct>mmand time* fhd polarity could be switched In phase 

with the motion cvery-othef dycle^ Switching everySther cycle makes 

the damping only 50^ as effective, but ltdVad still# bl^ Improvement 
• • • 
fvet the magn4ife hysteresis. 

• • 

We decided to Implement a digital phase- locked loop In the 
Alight computer to lock on to the phase of thj tumble motion and control 
the s-coll switching, there were some men-tvfvla^ pfoblems to overcome 
In this Implesientatlon.* First of al|m we were performing the computa- 
tlons on a computer that had fiaattng po|nC arithmetic capability 
(either hardware or software) and also had no hardware divide capabllltyt* 
As a consequence we had no programs fe# generating the trigonometric func- 
tions needed to censttuct /digitally | |kt local oscillator signal Inside 

the computer. To work around this difficulty we struck upon the Idea of 

• • 

locking a "sat^ tooth" functloit onto the magllet t ft p ftr signal, rather thqsi^ 

a sink wave. The unit amplitude saw- toothy 4lrdwn below can be generate^ 

• * • • 
fr#m ^simple fixed point operations. 




% 

%• 


••• 

e 




** 


Tha average delay Is based on #ne frame since the first frame Is 
etually likely to occur Just before or Just after the derivative 
changes sign. 

fits Hft^ Aamputer, 1 A keeping with Its primary functions, was oriented 
•toward tsgif and bit string manipulations rather than arithmetic 
operations. * 
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Locked on 



Figure 8. Saw Tooth Function for Phase-Locked Ijoop 


Ac any time t Che saw-tooth value, f(t) Is given by 
f(t) ^4tt)j(t • ^5% when co^ft - T^.) ^.5 

f(t) « 3 - 4u)j(t - T^), when uj^(t - T^.) i .5 

where 


(Uj " current saw4:9oth frequency • 1/T^ 


* tlBie of last saw-tooth trough. 

After estimating the amount of work Involved In d^el<^ng a 
software floating point package, we decided that the phase- locked loop 
could be written in fixed point arithmetic using the saw-tooth function 
The calculations had to be properly scaled, and double and triple word 
precision was used in some places where needed. 


The calculations were designed to avoid dlv^<«lon wherever 

f 

possible. In two^|||a<^s it was unavoidable, so we wrote a software 
division algorithm^ ing an iteration method for the inverse. This con# 
slats of Iterating the following equation: 
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vliere X ia the number «e with to obtain the Inverae of, and Y la the 
eat^aated inverse* The method converges for all positive X < 2, and ve 
insured convergence by proper scaling. 

The error signal for the phase* locked loop mas the product of 
the sampled elgnal (the a magnetometer reading) and the sampled VCO out- 
put (Hie internally generated saw-tooth function). The filtered error 
slgni^hen drove the saw tooth In phase, frequency, and frequency drift. 

The main problem to be overcome was the aliasing problem caused 
by the low sampling rate of the TM system. Our sampling frequency was 
.236 he and we were trying to stay locked to a signal (the -z magnetometer 
reading) whose frequency was In the range 0 to .75 hz. As the tumble 
rate decreased, the signal frequency passed through the imiltlples of 
the sampling frequency causing singularities where the error signal was 
driven to a constant value. We avoided this problem by making the loop 
of second order with variable gain, and using the frequency drift to 
"flywheel" through the singular points. As the frequency approached a 
multiple or half multiple of .236 hz, the gains were reduced to zero so 
that the saw tooth was running open loop. As t!te frequency passed 
through the singular point, the gains were gradually reestablished to 
their closed loop values. 

We^iere aided In the Implementation by knowing rather accurately 
the frequency drift when the z magMt was both off and on. These were 
determined experimentally early in the game using the DAMP program at low 
tumble rates. It can be shown from simple mechanics that the effect of 
the z-coll on the tumble frequency Is Independent of the frequency. We 
had determined the effect to be .09 rpm per 10 minutes of z-coll opera- 
tion near the earth's equator. Thus we could switch to the appropriate 
drift value as the z-coll was activated or turned off, and thereby avoid 
Introducing large transients into the loop. This meant we could duty- 
cycle the z-coll as necessary and still maintain phase lock. 

We could also determine to within about .1 rpm the frequency 
of the tumble motion by Independent analogue means. Since the tumble 
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notion rotatod th« apacacraft antanna, tha ground racalvar AGG voltagaa 
could ba nonltorad Co datamina tha Cunbla rata. Thla method was used 

to give tha program starting values for tha frequency when needed. 

a 

Va provided an AGC feature to conpensata for tha varying 
amplitude of tha z magnetometer signal caused by the change in the earth 
field strength with orbital notion. The field amplitude changes by a 
factor of 3 within one orbital revolution. This control was accomplished 
by filtering the absolute value of the signal and using it to normalize 
the magnetometer readings to unit amplitude before generating the error 
signal. The readings also had to be corrected for the effects of the z> 
coil dipole. 



The details of the loop are shown in Appendix B. We determined 
the appropriate gains and filter time constants by making a computer 
simulacion of the complete process and experimentally adjusting things 
until we could make it tfork through the range of frequencies from 
VS rpm to 2 rpm. found that the phase error would build up to 20° or 
so as the frequency passed through the aliasing points, but this was taore 
than enough accuracy to provide effective de- tumble. The simulation showed 
that the loop could achieve phase lock from a 180° initial phase error 
with an initial frequency error of .2 rpm. The time required to lock on 
was about 20 minutes. 


The system was Implemented in the flight computer as two separate 
programs: (1) a special version of the TM monitoring program \diich pro- 

cessed the magnetometer readings and locked onto the phase of the motion, 
and (2) a program operating off the clock interrupt which used the 
latest phase information from the first program to control the z-coil coji- 

t nds. The two programs ran independently, but worked in concert, passing 
fonsation back and forth through shared memory locations. The phase 
lock program ran continuously, and the switching program was cycled on 
for about 15 minutes each orbit period. 
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6. Gan«r«tlng « Tumbl* Notion for Centrifugal Force 

After the de>tuable wee complete, the gravity- gradient boom 
waa erected on HP- II* Becauee of an unforeaeen problem vlth aclaaora 
bocue, the link* broke during the deployment. It waa too lece to change 
the boom on the next apacerraft; but the problem was understood, and the 
boom could be successfully deplo]red If it wss kept in tension using' centrl- 
fugal force. This centrifugal force could be generated by tumblli^ the 
spacecraft at the correct rate. 

On HP-III the boom was successfull> deployed in this manner. 

The spacecraft was tumbled by using the phase- locked loon program and 
simply reversing the z-coll polarity commands. It wus important Jurlng 
the tumble-up operation that the rate not get too hl^ at any stage or 
the boom links could have been broken In tension. It proved easy using 
the existing special programs to let the flight computer automatically 
Increase the tumble momentum bit by bit as the boom was gradually 
deployed over a three-week period. 

V. FIHAL REMARKS 

As a result of the fll^t computer control capability, the 
TIP- II and TIP- III spacecraft achieved a substantial partial success. 

The TIP- III spacecraft was able to achieve 3-axls stability with the 
drag-compensation system In full operation. This allowed valuable 
in-orblt testing of this Important sub-system* The TIP-II spacecraft 
was able to go operational as a navigation satellite for limited periods 
(when the percent sunll^t was high enough). As the first spacecraft 
In a new series it also provided valuable training and debug capability 
to the operators of the Mavlgaclon Satellite Syatem. 

The experience proven to be a dramatic illustration of the power 
of having a programmable computer on board a spacecraft. By having enough 
flexibility built into the system, we had a powerful capability to modify 
the onboard logic after launch. We seemed to be limited only by our 
own Ingenuity In finding ways to use the system. 
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APPENDIX A 

FINDING AN OATS FIRING POINT GIVEN AN ATTITUDE 
We define the S, T, and W directions in Inertial space as 

follows: 

S in the direction of the radius vector to the spacecraft 
T normal to S in the orbit plane in the direction of notion 
W normal to orbit plane along angular momentum vector. 

The attitude of the spacecraft when spinning about its longi« 
tudlnal axis is given in terms of the ri^t ascention, or, and the declin- 
ation, 6, of the spin vector in the direction of thrust, a is measured 
from the equinox and 6 is measured northward from the equator. Hence, 
the unit vector in the thrust direction is 


^ /cos 6 cos cr\ 

F ■ I cos 6 sin crj 
\ sin 6 / 

^ A 

The STW components of r are related to F by 


(A.l) 



R, O) (i) R, (0) F 

Z X z 


where R (9) is a rotation of angle 6 about the "a" axis 

A 


(A. 2) 


0 is the argument of latitude 

1 is orbit inclination 

n is orbit ascending node 

He can define components A, B, and C as 

- R^ (1) R^ (n) F , (A. 3) 
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1** 

a* 


so that 


/ S \ / cos 3 sin 3 0 \ / A \ 

(Tj.(-.jn.c..0) (.) 


(A.4) 


Multiplying out the matrix equations A.l and A. 3 we get 
A “ cos 6 cos or cos Q -I- cos 6 sin or sin Cl 

B • -cos 6 cos or sin 0 cos i + cos 6 sin or cos H cos i -I- sin 6 sin i 
C • cos 6 cos or sin Cl sin 1 - cos 6 sin or cos 0 sin i 4- sin 6 cos i 

The planetary equations to zeroth order in eccentricity can be 
integrated assuming an impulsive thrust with direction cosiponents S, T, and 
W. Using Eq. A.4 the Integrated equations can be w: itten in terms of A, 

B, and C 


Aa ■ — I -A sln3 + B cos 3 1 FAt 


(A.5) 


FAt 

na 


A cos 3 sin f + B sin B sin £4-23 cos 3 cos f - 2A sin 3 cos f 
FAt 

Ai cos 3 “ C 

na 


where F > the magnitude of the thrust, 
f ■ the true anomaly 

These give the changes in the orbital elements, a, e, and 1 for a firing 
of duration. At. The thrust magnitude is known a priori from the fuel 
tank pressure, and the thrust duration is selected to be of significant 
length but not too long that the impulsive nature of the thrust is destroyed. 
Generally this is 4-8 minutes. 

The optimum changes in the orbit elements (Aa.^, 
calculated as if the same duration thrust were to be made, but with free- 
dom to chose the direction (a, 6) as well as the true anomaly. The equa- 
tions for the Aa^, Ae^ and Ai^ are not given here, but the results are 
such that the ratios 
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and 




are maintained constant from firing to firing. 

The program to determine the time to fire loops through the true 
anomaly In 1** increments completely around the orbit. At each position 
the changes Aa, Ae, and A1 are computed and a figure of merit constructed: 


F.M. 


/„ /4§. . , 

l\^ + w - 

if + „ . 1 ) 

1 iWi 




where the weights ^3 Inputtable. If the orbit changes 

were all equal to the optimum changes, the F.M. would be 0. The pro- 
gram determines the true anomaly which minimizes the F.M. , and repeats 
the process for a new set of Input weights on operator command. Thus 
the results could be obtained for various cases and compared before a 
choice was made. 


An example of the cases examined was to set W. ■ 0 and U, 


W, 


3 1 2 

so as to remove any constraint on changing the Inclination, and Just do the 
best we could at raising the orbit and circularizing. The Idea here was 
that If we could get a really favorable along- track thrust today at the 
expense of a small inclination change In the wrong direction, we could 
probably make up the Inclination change on the next thrust. 

Once the optimum true anomaly was selected, the program then 
calculated from the orbit geometry the length of time after set which 
would center the firing on that position. This was then converted to 
a "dummy" setting for the satellite clock, and displayed for the space- 
craft controller to use If the decision was made to fire. 


1 . 
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APPENDIX B 


DIGITAL PHASE- LOCKED LOOP FOR Z-COIL CONTROL 


Each TM frame, the program receives a sample of the raw z- 
magnetometer, , which we refer to as the ith measurement. This 
measurement is then processed to give an updated phase and period for 
the saw-tooth function to control the z-coil commands. The loop is 
shown schematically in Fig. B. 1. The computations are given in the fol- 
lowing steps. 

1. Each measurement, correct for effect of z-coil - 

m^ - - KA 

where m^ ■ corrected measurement 

“ raw measurement 

A * z-coll dipole moment 
K • jf 1, depending on z-coil state 

2. Normalize measurements - 

AMj^ - (1 - A) AM^_j^ + A|m^| 

""i 

FMj “ “ normalized measurement 

1 1 . !) AM 

A - 1/256 


3. Compute saw-tooth value based on last known trough time and frequency - 

■1 


where 


\ ' hi - V'] 


modulo 1 


(fractional part) 


jfj^ • 4rj^ - 1 , if r^^ < .5 
[fj^ ■ 3 - 4rj^ , if rj^ 2 .5 

r^ * fractional part of a cycle since last trough 
t^ * time of present measurement 
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4. 


5. 


6 . 


■ tixne of last trough 

uu ■ current saw-tooth frequency 
f^ ■> value of saw-tooth function 

Compute and filter error signal recursively - 

troduct) 

(Error signal) 


■ FM^ X fj^ (current product) 


*1 - ‘ ""^-l 


B 

C 

D 


3/256 

1/128 

31/32 


Check for aliasing points and reduce error signal if necessary - 

2 


If 


lU - U). 


a 


< a , - R-i 




6 ■ .1 rpm 

Apply controls to saw-tooth - 
tu « d) - K2 Rj^ 

0 ) • u) + (bAt - Kj^Rj 


(!iv^) 


"t - 'i ■ 


(U 


(drift) 

(frequency) 

(time of last trough) 


where r^ is the fractional part of a cycle calculated in Step 3, 

Commands to control the z-coil are then issued based on and 
l/«. The value of A) is changed to the appropriate theoretical value if 
the z-coll is turned off or on. This results in a small ramp type trans- 
ient. 
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A Kalman filter la |>reseate<| which processes pseu4h>« 
range and delta-pseudo-range data from the NAVSTAR Global . 
posiiioiiing System for the onboard navigation of a host satellite* 

The formulatton and structure oi prc^Bm components is described* . 
and sample results are given* 
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ORBIT DETERMINATION ACCURACIES 
« USING SATELLITE-TO-^ELLITE TRACigN^ 

F. O. Vonbun, P. D. Argentiero, P. E. Schmid 

^ Goddard SpaOe Flight Center 

1.0 INTRODUCTION 

The possibility of using geostationary satellites for communications was discussed in 
the popular literature as early as 1 956 (1). The first detailed proposal for a synchronoifs 
tracking satellite system for the purposes of orbit determination was provided by IfU^Aun ^ 
in 1967 (2, 3). Since then a number of papers (4, 5, 6, 7, 8, 9,) have considered the use of 
satellite-to-satellite tracking for orbit determination and for gravffy ffeld model refinement. * 
These papers mention that with regard to coverage, a sateUite*to-satellite tracking system 
has a significant advantage over ground based tracking systems. For instance, with a single 
synchronous relay satellite, a satellite-to-satellite tracking system is capable of observing an 
earth orbiting satellite during almost half of every orbit. Equivalent coverage of a satellite 
in a high inclination orbit would be difficult to obtain with a ground based system. 

In 1 968 during the early planning phases of the geostationary ATS-6 and near-Earth 
NIMBUS-S experiments it became clear that this satellite t^iinguratTon would be ideally 
suited to evaluate the concept of satellite-to-satellite tracking and to p^vide valuable ex- 

-i ' * 

perience in processing this new data type. The experiment as defiftpd in October 1968 (10) 
incorporated both radio time delay (range) and Doppler frequency shift (ranglme^niM*- « 
'foremtQts. Jhis experiment, entitled the ‘Tracking and Data Relay Experiment” (Td^lMU^* 

« „ e * 

° • o 

watoM^ctefoa planned oxcept that NlMBUS-6, which was launched June 12, 1975, 
rather fbaii NIMBtlS>$ milted fhe T&DRE equipment. In early 1972 plans were completed 

o • • ^ O 

0 , * 

for a very ATS-6/GE0S>'| eatellite-to-satellite tracking experiment. The GEOS-3 

O • * ' Q * 

O 4 * 

satellite wal launched on Apt9 9, 197S. Another satellite-to-satellite tracking effort involv- 

• A * “ ® ^ 

ingthv°AtS-6was the Goddard Af»flo*Soin>tCeodynamics Experiment (1 1) perfonned 
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during 1 97S. However the accent of this experiment was gravity anomaly detection rather 
then orbit determination. The ATS-6, which was the relay satellite for these experiments, 
was launched on May 30, 1974 and is still in operation. 

The results of these experiments are relevant because NASA intends to use the Track* 
ing and Data Relay Satellite system (1 DRSS) ( 1 2) for operational orbit determination of ^ * 
NASA satellites. The system will consist of two synchronous relay satellites (one at 41 
degrees west and one at 171 degrees west) and a common ground station under construction 
at White Sands, New Mexico. Operations will begin in November 1980. Hence by the early 
nineteen eighties sateIlite-to*satelIite tracking data will be routinely processed to obtain 
orbits. 

This paper is a report on the results of the ATS-6/GEOS-3 and the ArS-6/NlHlBUS-6 
satellite-to>satellite tracking orbit determination experiments. The tracking systems used in ^ 
these experiments differ from the TDRSS, primarily in the use of one rather than two syn- 
chronous relay satellites. However the authors believe and simulations mentioned in this 
paper indicate that the insights gained from the experiments with regard to proper data rcs 
duction techniques and expected results are applicable to the TDRSS. 
l.l EXPERIMENT SPACECRAFT 

• ^ ^ * The key to all satellite-to-satellite experiments to date has been the geostationary ** 

* 

AT|4^ spacecraft (13, 14). During the past three years the equatorial ATS-6 has Been sta- 

• 4 

tioned tiiih the Pacific in proximity of continental U.S.A. and Africa. 4»M:ordingly, 

» ATS-6 ^«fnd sfatiqns have at various times been operated at Rosman, North Carolina; 

Mojave, California 4o4^drid, Spain. The near-Earth satellites tracked via ATS-6 have been 

• * • • 

.* *. GEOS-3(IS),Apollo-Sdytti.M6|.^dNIMBUS-6(l7). 

• • ^ • 

The nominal GEOS-3.«rtni4l|!t»»ara«;|4^«iU'e adteatf altitude of 843 km. an inclina- 

* • • • * • ’ * • * * 

• • ' ***** •**•**• 

• « « * , tion of 1 1 , an eccentricity of 0.004:*4il«t J Ntiod of lOLt mn, IW 

' - ' • • . 
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were chosen to minimize resonance of the subsatellite trace with any given Ear^tfedhire 
and to provide orbit traces which fbver the Earth in a gridwork pattern. 

The Apollo-Soyuz mission included the Geodynamics experiment where Apollo was 
tracked via ATS-6 for the mission duration (15 July to 24 July 19?S). The nominj||ppMlo 


orbit at insertion was ISO km by 170 km at an inclination of SI. 8°. 


♦ ® «v 
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Finally, the NlMBUS-6 weather satellite is in a Sun synchronous polar orbit with a« * « 
mean altitude of 1 1 10 km. an inclination of 100 , and a period of 107.4 min. i 

2.0 SATELLITE-TO-SATELLITE TRACKING 




Sfatellite radio or laser tAsking system makes measurements of such parameters as 


range, range rate, angles and direction cosines to a spacecraft relative to a giv#N tracking sta- 

9 

tion. In two-way tracking a signal i$ transmitted from a well surveyed ground station to a 

spacecraft transponder which frequency translates the signal for re-transmission , 

* 

back to the ground station or, as in the case of satellite-to-satellite tracking, to another 

spacecraft. The two-way tracking system developed for the experiments discussed in this 

« • 

paper measures “range” in Icfms of the round-trip tmve delay on a 100 kHz tone and range 
rate in terms of the Doppler shift oh « 2 GHz carrier signal (14, 18). 

2.1 GEOMETRY ^ 

The tracking geometry is shown in figure i. Tlie ground sta a transmits a signal to 
the near Earth satejlitc via the synchronous spacecraft. This same signal is “turned around" 
and transmitted (at a sliglitly offset frequency) back to the ground site again via the high 
altitude satellite. For purposes of stability NASA geostationary orbits have been maintained 
at inclinations which extend from 1.5° to 6°. As a consequence the path indicated as Rj 
‘ (figure i) varies a§a function of time as docs R^ (13). Because of the radio propagation 

• -.times involved and the fact that both spacecraft are in motiem relative to the ground .site, 

• four distinct paths must he considered when interpreting the Doppler (range-rate) and lime 
' del'ayrange/measufements<-|3, 19). 


V « 




SYSTLM DESCRIPTION 

•• • 

“range" measurement is performed by comparing transmitted and received tone 
zero crossings, the highest resolution tone frequency in^liis case being 100 KHz. Lower fre- 
quency tones are sequer ially used during acquisition for ambiguity resolution. The lower 
tones are at 20 KHz. *4t4iz^«CO Hz. IdCJU i2 Hz. and <S Hz. 

* The tone ranging .'*%';i*urement is qK^e4Jralghtforward and ranging accuracy depends 
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chietly on the quality <of prellight calibration of both^n A^'S and Nl.MBUS transponder 
o o ®o • 

group delay. Such v^elTiglU caliRraJfeite.irfae'a have been lake^t over a riige of frequencies and 

o ^ ® • 

iemperatures. lndicat>i«n^ate®t#f !<«riffee'»-i»^eful calibration the%i^ksystematic delay error 

in the '.anging measureinein ci^fbe hel#4;oa hew meters of eqftvilent one-way range. 

6 ® * 

On • 

The “range rale" measure'ncnl^4i®etitffmed by eountiig cycles of Doppler over a 
measured time interyl. In the case l.OS-3 ami Apollo the measure!nent consisted of 


the number of Doppler cycles accumulated in the regular sampling inlerval ( 1 or 10 seconds 



depending on mode of operation). For N1MBUS*6 the measurement consisted of the time 
interval required to accumulate a fixed number of Doppler cycles (18). The electronics 
for ATS^ satellite-to-satellite tracking have been so configured that the Doppler output is 
approximated by: 

-2f,k _ _ _ 

fd=-~Uir*,+a2(f,+f2)l 

where 

fjj = measured average Doppler frequency 
f| = uplink frequency 
c ® speed of light 

k, a I and a 2 are scalar constants determined by equipment frequency multiplications 
#, = average range-rate ATS to ground site 
?2 = average range-rate ATS to NIMBUS, Apollo, or GEOS-3 
A detailed discussion of Doppler factors in satellite-to-satellite tracking is given in (19). 

The uplink to ATS-6 (f^) is at a nominal 6 GHz. The link to and from the low satellite is 
nominally 2 GHz and ATS-6 back to ground at 4 GHz. The range and Doppler measure- 
(P ments will also be biased by the Earth’s troposphere and ionosphere. Measurement biases 
^p to meters in range and tens of cm/sec in range rate can be expected at 2 GHz. Atmos- 
phere refraction effects can to a large extent be modeled out. Some of the work done in 
this area at NASA-GSFC is indicated in (20, 2 1 , 22). The atmospheric range bias is fre- 
quency independent through the troposphere and inversely proportional to frequency 
squared through the ionosphere. The range rate bias, in addition to the foregoing, is pro- 
portional to the rate of scan through the atmosphere as well as to the magnitude of hori- 
zontal gradients. 
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2.3 ORBIT DETERMINATION TECHNIQUES 
Problems of Orbit Determination with Satellite>to^atellite Tracicing 

Tlie unfamiliar feature of determining user satellite orbits by means of a satdlite>to- 
satellite tracking system is the presence of the relay satellite state as an error source. The 
simplest procedure for estimating user satellite state in the presence of this error source is to 
estimate a satellite epoch state from the satellitC'to>satellite tracking data with the relay 
satellite state constrained to a previoudy determined orbit and left unadjusted in the reduc- 
tio'i process. With this approach the uncertainty in relay satellite state is manifested as an 
unmodeled and time varying error source which alters the estimate of user satellite state. 
Some .subtleties are encountered in determining the effect of this error source. The time 
history of relay satellite state error is a function of the way in which the epoch state was 
computed. For example, suppose the relay satellite is independently and continuously 
tracked over a given period and a least squares algorithm used to estimate epoch state at the 
I'ecinning of the period. If this epoch state is then propagated to the eiid of the period 
using the same dynamic model that was used to process the data, the resultant errors will be 
constrained by the data fitting criterion implicit in the least squares reduction algorithm. 
The errors so obtained will be smaller than the errors obtained if either one did not match 
'.lynamic models or if one propagated the epoch state beyond the data collection period. 

The same phenomenon can be understood from a statistical vantage point by observing that 
when the dynamic models are matched the epoch state errors become correlated with dy- 
namic parameter errors, and that over the data arc these correlations tend to minimize the 
errors in the epoch state propagation. 

Relay satellite state uncertainty appears to be a significant error source even when the 
relay satellite or satellites are continuously and independently tracked. Argentiero and 
Loveless (23) simulated the orbit recovery of a satellite in a 300 km, polar, circular orbit 



with the Tracking Data Relay Satellite System (TDRSS) (24). The TDRSS satellites can 
relay range and Doppler information from a low altitude user satellite to a ground station. 
The simulations assumed that each synchronous satellite was continuously tracked from 
two ground stations and that 24 hour data spans were processed to estimate user satelUte 
state. The same dynamic models which were employed to estimate relay satellite epoch 
states were also used to estimate user satellite state from the satellite-to-satellite tracking 
data. The effect of Geopotential and atmospheric drag errors were included in the simula- 
tion. The results showed that user satellite position could be recovered with an average 
total position error of 260 m. The major part of this error is caused by the error in esti- 
mates of relay satellite epoch states. When these simulations are repeated without the as- 
sumption of continuous tracking the results are considerably worse. 

A standard approach to dealing with troublesome error sources in an orbit determina- 
tion is to augment the list of estimated parameters in the data reduction by including these 
error sources. This approach can certainly be implemented with regard to relay satellite 
state errors by simultaneously estimating user and relay satellite epoch states from informa- 
tion supplied by the satellite-to-satellite tracking data. From one vantage point this is an 
undesirable solution in that the user is uninterested in the state of the relay satellite and 
would rather not burden the numerical procedures with the need for simultaneously esti- 
mating relay satellite state with the user satellite state. However, the results of independent 
covariance analyses performed by Fang and Gibbs (2S), and Ai^entiero and Garra-Robles 
(26) indicate that an unconstrained simultaneous estimate of user and relay satellite states 
using satelUte-to-satellite tracking data can yield an estimate of user satellite state which is 
consistently better than 100 m. 
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Numerous simultaneous unconstrained solutions have been attempted using range 
sum and range sum rate measurements obtained from the ATS-6/GEOS-3 combination and 
the ATS^/NIMBUS-6 combination and in all cases the solutions have been inaccurate and 
numerically unstable. Clearly our experience with real reductions of satellite-to-satellite 
tracking data is not compatible with the results of previous error studies. In order to under- 
stand the discrepancy we have performed a numerical simulation of the ATS-6/GEOS-3 
satellite-to-satellite tracking experiment. The difference between a numerical simulation 
and a covariance analysis can be described as follows: in a simulation, data are generated 
and a least squares adjustment process is actually performed. The estimated state is then 
compared to the reference or unperturbed state at various points along the orbit and conclu- 
sions can be drawn concerning the accuracy of the process. In a covariance analysis mode, 
the least squares adju.«tment process is postulated rather than actually performed, and under 
the assumption that over the range of expected errors, perturbations of orbital estimates are 
approximately linear functions of perturbations of the error sources, the associated covari- 
ance matrix of the epoch state recovery is computed. With the aid of state transition matri- 
ces the covariance matrix at epoch can be propagated to obtain the covariance matrix of the 
satellite state recovery at any point in the orbit. 

For the numerical simulation a computer program was used to generate 1 2 hours of 
range and Doppler satellite-to-satellite tracking data from the ATS-6/GEO.S-3 satellite com- 
bination. In this data generation the Naval Weapons Laboratory (NWL) g«opotentiaI field 
was used. A random number generator added white noise of standard deviation 1 mm/sec 
to the Doppler data and white noise of 2 m to the ranging data, values consistent with track- 
ing system performance. The SAO 69 geopotential field and an orbit determination pro- 
gram were used to reduce the data to simultaneously cstim ; the ATS-6 and GEOS-3 epoch 
states. The estimates GEOS-3 epoch state was propagated along the entire 1 2 hour data 
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collection period using the SAO 69 geopotential field. This orbit was compared at selected 
time points to the true GEOS-3 orbit which was obtained by propagating the GEOS-3 ref- 
erence epoch state with the NWL geopotential field. The average difference between the 
two orbits was over 900 m. Also the nominal covariance matrix of the data reduction re- 
vealed that several correlations between estimated parameters were of absolute value near 
unity. This implies that the normal matrix which is inverted in the least squares estimation 
process is poorly conditioned. Hence small perturbations of the elements of this matrix 
such as those caused by computer roundoff and other effects cause m<uor perturbations 
of the elements in the inverted matrix. This amplification effect in the inversion of a poorly 
conditioned matrix can lead to an inaccurate estimate of a satellite epoch state or in some 
cases a divergence of the least squares interation procedure. This is the probable cause of 
poor results using a simultaneous estimation approach in both the simulated and real data 
reductions. In a covariance analysis of the simultaneous estimation approach the least 
squares algorithm is not actually executed and consequently these numerical problems are 
never manifested. For this reason the techniques of covariance analysis provide a somewhat 
optimistic assessment of orbital accuracies obtainable from simultaneous estimation with 
satellite-to-satellite tracking data. 

Thus the two conclusions of our analyses are: 1) The uncertainty in relay satellite 
state is a significant error source which cannot be ignored in the reduction of satellite-to- 
satellite tracking data and 2) that based on both simulation^ and real data reductions it is 
numerically impractical to use simultaneous unconstrained solutions to determine both 
relay satellite and and user satellite epoch states. The estimation technique used to generate 
the results shown in subsequent sections may be described as a Bayesian or least squares 
with a-priori procedure. This approach permits the adjustment of relay satellite epoch state 
in the reduction of satellite-to-satellite tracking data but without the numerical difficulties 
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introduced by an ill-conditioned normal matrix. Theoretically this technique obtains the 
best possible estimate of user satellite state based on all available information. A mathemati- 
cal description follows. 


Mathematical Description 

In this mathematical development we assume the existence of two separate data sets: 
y I - Ranging observations between ATS-6 and ground based tracking stations 
y-) - Satellite-to-satellite tracking of user satellite (range sum and range rate sum) with 
ATS-6 as relay satellite. 

The parameter set to be estimated consists of two satellite epoch states. 

X I Six dimensional ATS-6 state at epoch time T| 

Xt Six dimensional user satellite state at epoch time Tt 

The data set y j is corrupted by errors in the measuring process. Hence represent yj as: 

V| =Vi +Vj.e(Vj) = 0.e(v,Vj^) = g, (I) 

where y j is the correct or noiseless representation of the data set, vj is a vector of random 
errors of zero expectation and covariance matrix Q| . Describe the functional relationship 
between y I and X| as 

y,=f(X|) (2) 

The right side of eq. 2 represents a computational algorithm obtained by integrating satel- 
lite motion to each observation time and computing the ideal observations. The standard 
least squares estimator X| of X| is that vectoi which minimizes the loss function. 

L(xi) = (yi - f(Xi))^Qi"'(yi - f(x,)) (3) 

Assuming the linearity of eq. 2. tlie vector which minimizes the right side of eq. 3 is also 
known to be a minimum variance estimator. A first approximation to the desired minimum 
can be obtained by expanding equation 2 in a first order Taylor series about nominal value 

’‘!>n 
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( 4 ) 


where 6yj and 5X| 
sensitivity matrix. 


6yj = Ajfixj, Aj = 


9f(xi) 

3x. 




are deviations of 7] and Xj from nominal values and Aj is the so-called 
The estimate of 6xj is 


6x, =(A,TQ,-lA,rlA,'^Q,-*5y, 


( 5 ) 


where 

«yi=yi-f(xi,„) 

The vector 8x j is added to x j „ to form an estimate of x j . This estimate can be used as 
a new nominal and the process can be repeated until a convergence criterion is satisfied. 
The covariance matrix of the least squares estimate Xj of Xj is 


c = e((x, -x,J (Xj -x,]T) = (ATQ,-lAr‘ 


( 6 ) 


The next step is to obtain an optimal processing of the data set Define a 12 dimensional 
vector z as 


z 



( 7 ) 


Represent the data set 


y2 = ?2 + ^2* * Q 2 


( 8 ) 


where y 2 is the correct or noiseless representation of the data set, V 2 is a vector of random 
errors of zero expectation and covariance matrix Q 2 . The functional representation be- 
tween y 2 and z is presented as 

y2 “ 8(z) 
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As was t!ie case with equation 2, the right side of cq. 9 a'presents a computadkal algoritiun 


involving the integration of satellite equations of inotio“ 

0 

The least squares estimate of z would not be optimal unless all available information were 
included in the loss function. Hence it is appropriate to treat the least squares or minimum 
variance estimate X| of Xj as an a-priori estimate weiglited by the inverse of the covariance 
matrix provided by equation 6. The resulting toss function to be minimized is 


L(z) = (y-, - g(z))^0-,"’(yi - g(z)) 




Again, the required minimum can be obtained iteratively by expanding equation eq. 9 in a 
first order Taylor series about a nominal estimate z„ of Z 


Sy, = A,6z. 


3 z I z=Zn 


where Sy-, and 8 / are deviations of y, and z from nominal values and A-, is the sensitivity 


matrix. The estimate of 5z is 


6/. = (aJq,-'^+c"' 


A-,^Q2“*6y, ^ ® “ ] 

1^0 c~y LSx.J / 


where 


5x, = X, - Xj , 5y, = y, - g(z„) 


The vector 6z is added to estimate z. This estimate is used as a new nominal and the 
process is repeated until a convergence criterion is satisfied. The final covariance matrix for 
the estimate of satellite state Xj and satellite state x, is 


= { + c“^ 


V 


• •• 


It can be shown that the two step process defined above in which data set yj is processed 
and then data set y 2 is processed is equivalent to a single step unconstrained least squares 
estimation of Xj and X 2 using both data sets yj and y 2 < Hence this procedure leads to the 
most accurate estimate of both user and relay satellite state based on available information. 

3.0 EXPERIMENTAL RESULTS ^ 

u 

The experimental results can be consider'd in three categories - namely, 

A • 

• tracking system performance ^ ^ 

® A 

• geostationary satellite orbit evaluation ^ ^ 

• near Earth satellite orbit evaluation 


o o 
• o 
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The expected error for the NAS 4 range and range rate satellite-to*satelliteftracking 


3.1 TRACKING SYSTE^ERFORMANCE 






9 

^stem is a function of many controlled parameters^uch as range tone frequency; S^ple 

^ ® ft ® * 

rate, bandwratn settings, signal<to>noise smctral density ratios, spacecraft dyn^iCis,^^ O 


o® • 


P P 

O ° . 

Q Op# 

P 

o 


(13). However, the system is generally used with what |iight be termed a standard set of 

A # p p 

options sucmS: 100 kHz maximum range tone frequency|pignal levels such that system is 
not thermal noise limited, 1 per second or 6 per minute data rate, and a 2S Hz range track* 

O • 

ing loop two-sided noise bandwidth. Table I lists the theoretical system performance for 
the foregoing selected options. Doppler averaging tin^s appr^imately one h^ the sample 
time interval for NIMBUS tracking and equal to the sample intemQfor Apollo 4 u^E 0 S 

P 

tracking. H 

For averaging times, T, up to about 10 seconds the noise d^eases as^T. Che prin* 

P • ^ • • 

cipal Doppler noise contribution comes from receiver voltage controlled crystal •scillators 

• * • 

and the analog to digital conversion. For longer integration times the Doppler noise is also 


ft Oft 


influenced by noise falling off as 1/T, an effect attributed to t|e phase jitter in the 
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transmfttep reference signal use^ at the Dopplft* extractor. For safellita-to-sa!fllileIfj*cking 

^ ^ « u‘ • '- 

.r> o 

itivolving ATS-6 the range-rate resolution for T seconds of averagin^(27^ is given €>y: 
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This range-rate j?^solution veriiv yopploJE measurement averaging time is «>Figure 2. 
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^It should be mentioned that tBo|?ast sigi)|ficant range bit recorded is %,5 meters 
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whiebis consisfcnt witb the besfcexjJtcted one way perform ince of 1.7 m»to*» rc«qJution. 

^ # 

asured results indicate close a|reenient with exported system performance. System 
random errors or “noise’'4irc geimrally observed hi' the least squares fitting ©f slfcort flata iP 

ill ^ 

spanifi.e.,^ tb 10 minutes) with polynomials of Jf least 5th degree to account f«r^aee* 

craft dynamics. Care must be taken such that the polynomial itself does not introduce 

s 

apparent error. If the data is from a static or colMmatfcn tefcver test a least squm’s straight 
line fit is appropriate. 

Assumitig reasonable tracking geometry the accuracy of spacecraft position and ^ 
velocity determination will be primarily limited by tracking system,f»erformance for any 
computation spanning the data collectioi interval. That is. if continuoini trackit^ is pro- 
v^ed from a set of well surveyed statitrlf the computation is essentially one of getmetry. 
On the other hand the accuracy of orbit prediction based on an initial space, laft vector 
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to deittawM Itie t%pecit4 itnon. in the cilunatc of ATSh 6 «ialc frmn gfoond ttaned IracUtig. 
Mofe.f>r€<i»ely, we require antwen to {theie questions; , o ' 

1^ Mow accurately can tfcr ATS-6 orbit be r!etennine4 over an oi^lat pcfiod . .. 

• * * • • < 

(24 hr) from data which spans the orbital period? 

* ' 2) Once an ATS-6 epoch state is determined, how accurately lUp that state ^ ^ 

. * 
propagated beyond the data arc whidh was need in its estimation? 

GEOSTATIONARY SATELLITE SHORT TERli ACCURACY 
Ttur-first question was investigated by examining reductions of ATS-6 tdlateration 
tracking obtained on November 3, t974 Tritaleralion data iisobtaiiKd bf igttJiwe a signal 
from a single tracking station to several strategically deployed unmanned low cost trans- 
pmider* via' the satellite whose state is to be detennined. The time required tor the radio 

aignals to complete the round trip to and from each tramppnder is measured at the trans- 

• « * 

milter site-. Thejnterrogating sites were located at Rosman, North Carolina and Mojave, Cali- 
fornia. The transponders w*‘re located at Rosman, Mtijave, Greenbelt, Maryland, and San- 
tiago, 'Chile. * ‘ * 

« • 

The method of “orbit overlaps** war died loevj^hule tlr orbit detennination accu- 

racy of the system. This procedsiro can be outlined aa follows; 

• • • • . * 

. I| Octyrmhic a satellite epodi state using eacbof two independent data sets 
• « • 

. • • * * 

*'* ’ 2) Propagate estimated epoch states oyer a esnn^son or overlapping interval 

C 

' 'Tf Differcnee the two orbits over tM'omnitom interval (differences are usually 

R - * * 

* displayed in akmg track, ereifer Ira^k, and fndial components). 

O 

In nouH' case* the qrbit overlap mctiiod can lead to an undcr-cstimation of orbit errors 
Hney biases in snbii e<Uimates may cftivtl in orbit differences. Hence, the method should be 
vts'wrd as a teM of the internal consistency of an orhif determination process rather than an 
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absolute measure of Accuracy. Data set ! used in tte oH4t overlap test was obtained with 
Jtosman as a transmitting site and with transp>jnders at Rosman, Mcjave, Greenbolt, and 
Sainbago. Data set 2 was obtained with the same transponder sites but with the transmitter 
• located at Mojave. The tracking schedule is shown in figure 3. The two interrogating sites 

<9 ' 

*« 4 o . • 

are identified in figure 3 under TRANSMITTER as Rosman, North Carolina and Jhe Mojave, 
* California “Hybrid Transportable” station. Each data stretch was approxknateIy-5 minutes 

long and the data rate was one sample per 10 seconds. Separate orbit arcs were computed 
from data set 1 and data set 2. The total position differences between the two orf^tt over 
the 24 hours of Nov 3, 1974 were computed and are displayed in figure 4. The mean posi- 
tion error is about 100 m. A typical set of range residuals is shown in figure 5. The range 
residuals over tliis arc are on the order of 20 m. ’ 



UN'VERSAL TIME (HOURS) 


Figure 3. Tracking Schedule 3 November 1974 

Assuming that there are no significant biases in the trilateration orbit determination 
whose effects cancel in the orbit overlap test, the results of figure 4 suggest that continuous 
tracking of ATS-6 over a 24 houf period leads to an Orbit (Hrtimate over the period which is 
accurate to about 100 m. 
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Figure 4. ATS-6 Total Position Error 
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Figure 5. ATS-6 Orbit Residuals 
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GEOSTATIONARY SATELLITE LONG-TERM ACCURACY 

In general one cannot assume that relay satellites are continuously tracked. Hence, in 
the reduction of satellite-to-satellite tracking data, it may be necessary to use an estimate of 
user satellite state obtained through a propagation that was unconstrained by the data fitting 
criterion of a least squares algorithm. When this occurs the accuracy of the orbit estimate is 
entirely dependent on the correctness of the force models used in the propagation. 

The orbit overlap technique (30), utilizing data obtained during July 1975 was used 
to estimate the accuracy of a free or unconstrained propagation of an ATS-6 epoch state. 
The data sets used in the overlap tests weie; 

DataSet I — 24 hours of data over July 13,14,1975. fracking f.iuiions located at 
Madrid, Ascension Island, and Johannesburg. 

Data Set 2-24 hours of ranging data over July 25, 1975. Tracking stations located 
at Mi.Jrid, Ascension Island, and Johannesburg. 

Each data set was processed to estimate an ATS-6 state vector for epoch time July 
16, 1975 at 7 hr., 25 min. The epoch states were propagated forward for 10 days and along 
track, cross track, and radial differences were computed at 15 minute intervals. The root 
mean square along track difference was over 2 km. Figure 6 is a plot of these along tn^k 
differences. 

The large errors which occur during the free propagation of an ATS-6 epoch vector 
must be caused by a misrepresentation of force models. The obvious candidates are: 

1 ) Unmodeled venting and thrusting of ATS-6 to accomplish satellite attitude cor- 
rections. Motions due to antenna maneuvering may also introduce errors. 




4.5 I I I I I I I I I L. 

7 33 59 85 111 137 163 187 214 240 

HOURS 


Figure 6. Along Track Orbit Differences for ATS-6 Overlap Test, July, 1975 

2) Mismodeling of solar radiation pressure. In all data reductions. ATS-6 was as- 
sumed to present a constant cross section to the sun. In fact, this is not the case. 

3) An error in representation of the central force term of the Earth’s gravity field. 
An estimate of the uncertainty in estimates of this parameter is one part in 10^. 

Error source number 3 appeared to us as the most likely cause for the major part of 
the errors exhibited in figure 6. In order to measure the effect of uncertainty in the gravity 
field parameter on the free propagation of ATS-6, the following simulation was performed; 
ranging observations to the ATS-6 from sites at Rosman, Santiago, and Mojave were gener- 
ated for a three day span. The observations were corrupted with white noise with a standard 
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deviation of 10 m. The value of the gravity field parameter used to generate the data was 
perturbed by one part in 10^ and this value was used along with a least squares estimator to 
estimate an epoch state at the beginning of the three day data span. The perturbed value of 
the gravity field parameter was used to propagate this epoch state for six days. Over the 
three days covered by data the propagated orbit differed from the assumed true orbit by 
about 100 m. But at the end of the six day propagation period the errors were approxi- 
mately 2 km. The results of this simulation suggest that the error in the central force term 
of the Earth’s gravity field is sufficient to account for the errors in the ATS-6 free propaga- 
tion as manifested in figure 6. 

SUMMARY OF RESULTS 

Overlap tests perfonned with real data together with simulation results suggest that 
by processing data over one ATS-6 orbital period, the ATSr6 state over the orbital period 
call be determined with an average accuracy of about 100 m. But other results show that 
when longer data arcs are used or when an estimated ATS-6 epoch state is propagated well 
beyond the data arc used in its estimation, errors in the kilometer region are encountered. 
These facts indicate that there are significant errors in the models of the forces acting on the 
ATS-6. The most likely candidate is the error in representation of the central force term of 
the gravity field. 

3.3 GEOS-3 ORBIT DETERMINATION RESULTS 

The GEOS-3 orbit determination results were derived from data obtained over the 
weekend of May 3, 1975. The tracking schedules and the tracking systems used in the eval- 
uation are shown in figure 7. The figure shows that five passes of range sum and range sum 
rate data were available. A Bayesian estimation technique described in a previous section 
was used to obtain two separate and overlapping GEOS-3 orbits. A GEOS-3 epoch state at 
May 2, 22 hr was estimated using all the ATS-6 ranging data and the first three passes of 
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Figure 7. Data Summary 

range sum and range sum rate data. The ATS-6 ranging data was weighted aeeording to a 
standard deviation of" 2 m. The range sum and range sum rate data were weighted according 
to standard deviations respectively of 2 in and 1 mm/sec. The complete GEM-7 geopoten- 
tial field was used in this and all other data reductions. The estimated epoch state was 
propagated to the end of the data span of May .T 22 hr. The process was repeated with the 
last four passes of range sum and range sum rate data to estimate a CiF.OS-3 epoch state at 
Ma\ 3. 10 hr. This epoch state was propagated to the end of its da'a span at May 4, 10 hr. 
The total position difference between the two orbits during the 1 2 hr overlap period as 
shown in figure 8 varies periodically between 10 and 25 meters. 

As mentioned in a previous section, orbit overlap procedures can provide an overly 
optimistic assessment of orbit determination accuracy. A more objective measure of accu- 
racy is obtained by comparison with an orbit derived from an independent and well 
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calibrated data set. Figure 7 displays the C-band tracking available from Wallops Island and 
Bemuida during the weekend of May 3, 1975. A three-day (IFOS-3 arc was derived from 
the C'-band data and compared to a similar arc derived from the five passes of sudellite-to- 
satellite tracking data and ATS-P ranging data. The root mean square differences in the two 
arcs were; 

radial 5 m 

cross track 200 m 

along tiack 39 m 

Various orbit results indicate that total position error for C-band derived GEOS-C orbits is 
on the order of 50 m (31). Hence, it is only in the cro.ss track direction that the orbit deter- 
mination derived from satellite-to-satellite tracking data differs significantly from the 
C-band orbit. The large cross track errors can be explained in terms of the tracking geome- 
try. For each of the five satellite-to-satellite tracking passes shown on figure 7, the GEOi-3 
satellite passed ahuost directly under the ATS-6 satellite. Consequently there was little 
cross track information in the rainge sum and range sum rate data. It is a reasonable assump- 
tion tiiat with a better geometric distribution of passes the cross track errors \x’ould be sub- 
stantially reduced. 

3.4 NIMBlJS-6 ORBIT DETERMINATION RESULTS 

The NIMBUS-6 overlap results were derived from data obtained over the weekend of 
Feb 8, 1976. For this experiment a higlily accurate reference orbit suitable for the puirtose 
of comparison was unavailable. This implied that the primary measure of the quality of the 
orbits derived from satellite-to-satellite tracking would be obtained from orbit overlap test. 
Hence the orbit overlap test for the ATS-6/NIMBUS-6 experiment was performed in a way 
which was more rigorous and less optimistic than the overlap test performed for the ATS-6/ 
GEOS-3 experiment. Notice that for the ATS-6/GEOS-3 experiment the overlap test was 
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porl'ornied with data sets which intersected through the entire overlap interval. Hence both 
orbits used in the comparison were constrained by data at each end of the interval. With 
such a procedure it is possible for the effects of errors in the measuring system to cancel in 
the test results. It will be seen that for the ATS-6/N1MBUS-6 overlap test the two data sets 
in question are abiding rather than overlaping and effect of measurement system errors are 
less likely to cancel in the test results. 

The tracking schedules and the tracking systems used in the evaluation are shown in 
figine 9. The first two rows of this figure show the tracking schedules for the ranging data 
from Madrid, Spain to ATS-6, and from Ahmedabad, India to ATS-6. The third row 
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GMT HOURS OF 3 MAY 1975 

Figure 8. (iFOS-3 Overlap Position Differences 
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RECOVERY OF SECOND NIMBUS EPOCH VECTOR. 

Figuie 9. Measurement Periods and Period for Nimbus Orbit Comparison 

indicates the tracking schedule for trilateration data (Madrid-ATS-6-Ascension Island). The 
fourth row shows the tracking schedule for the range sum and range sum rate data with 
Madrid as the ground station. The vertical bar located at 9HR/UT Feb. 8 indicates the 
epoch time for two estimated NIMBUS-6 epoch states. Epoch state 1 was obtained by exe- 
cuting a Bayesian least squares estimator with all the ATS-6 ranging and trilateration data 
and all the satellite-to-satellite tracking data to the left of the epoch time. Epoch state 2 
was obtained by repeating the procedure with the satellite-to-satellite tracking data to the 
right of the epoch time replacing the data to the left of the epoch time. The horizontal bar 
in row 5 of the figure displays the common interval over which the two NlMBUS-6 epoch 
states were propagated. The complete GEM-7 gravity field model was used in all data reduc- 
tions and propagations. The relative weights for the data types were obtained by first using 
nominal weights and processing all the data to estimate ATS-6 and NIMBUS-6 epoch states. 


145 










Tlie residuals of the estimation were used to determine the standard deviations of the noise 
on the various data types. These standard deviations were used to obtain weignts for the 
final data reductions. The computed standard deviations are shown on table 2. 


TABLE 2 

Standard Deviation Used for Weighting Measurements in Nimous-6 
Satellite-to-Satellite 1 racking Orbit Determination 


RANGE (INDIA) TO ATS 
RANG^HdRID) to ATS 

168 m 

50 m 

TRILATERATION (MADRID. ATS. ASCENSION) 

15 m 

SATELLITE-TO-SATELLITE RANGE 

11 m 

sati;lliti:-to-satelliti; range rate 

.3 cin/sec 


The data reductions were complicated by the fact that an experiment onboard the 
NIMBUS-6 was responsible for some outgasing. This effect was modeled as a constant along 
track thrust whose magnitude was estimated in the data reductions. The recovered magni- 
tude was approximately 10“^ m/sec-. 

F'igures 1 0. 1 1 . and 1 2 display the along track , cross track, and radial differences in 
the two epoch state propagations during the overlaping period. The R.M.S. differences are 
40 m along track. 30 m aross track, and 1 2 m radial. The secular growth of residuals in the 
along track direction is explainable in terms of gravity field error and an imperfect modeling 
of the outgasing effect whose direction was probably not exactly along track and whose 
magnitude was probably not constant. 

Finally it should be mentioned that a NIMBUS-6 orbit derived from satellite-to-satel- 
lite tracking data was compared to a NIMBUS-ti orbit derived from minitrack data. Tlie 
orbit differences were well within the stated accuracy for minitrack orbits of 500 m. 
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4.0 CONCLUSIONS 

The ATS-6/NIMBUS-6 and ATS-6/GEOS-3 satellite-to-satellite radio tracking syste 
performs as specified with a resolution of 1 meter in range and .03 cm/sec in range-rate for 
a 10 second averaging. 

A Bayesian least squares estimation technique utilizing a good a priori estimate of 
relay satellite state was used during these experiments to obtain user satellite orbits with 
accuracies comparable to what is obtainable from ground tracking systems. The limituig 
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HOURS OF FEB 8, 1976 

Figure 1 1. Cross Track Differences for Nimbus-6 Satellite-to-Satellite 
Tracking Orbit Determination 

factor in an orbit determination with satellite-to-satellite radio tracking appears to be the 
accuracy of the force models rather than tracking system precision. 

The results of these experiments imply that with the proper data reduction proce- 
dures, the tracking data relay satellite system should provide orbit determination capability 
comparable to what is now obtainable from ground based systems. 


# 


148 



• •- 




• •a 1 ^ ^ ^ 




• •• . • 

• • ^ 


• ft 


• •*«.* 


•%. • 




• 8 


•is 

Ul 

9 *.-4 

• 9 


g OH •# • . 


• !• 




• ft • 


•ft -!• 


• • #«» «»•. 


^®3 Sft .« ft ’•. ft« •• 


• 10 *,$1 ,12 13 H* 18 


. . •* 

.ft ® 


„• • •• * woufis (fp P 6 B 8 . 1976 • • 

•• 


Figure 12. Radial ltifferences#er Nimbus-^ SateII||e*to>Satellite ’ 

® o * ® fracking Orbit P^terminatioai * ^ • • * 




• • • 

*• • #• 

ACKI10WLEt>0MEl«fS^ 


• !* ft 


" .**• • .• 


• • •• • 


• • «» 


- to • 


• - O 


O • 


^ 0 • o •• * ^ • e? . ft 

^|The author* lhank Cfr.Ge<Jrg Mopiuch*and Mr. fblim^ryaj? for thetf •*. •* ^ * *, • 

• * •* • *» * •*,.*•• •*. 

course of this«e%arch. • • **•.**•*.*.• 






• ^ .ft • . 

• • • • 

• • • • • 

• « . • 

. • 


»«r 



'■ .. REFE«fcNm .• 

• li *‘^ 'Serf«tdl*K,iiidW. Boiler, “Srfrlli*)?”, f*. IJ5, tl4nov^Houi<.€*srd<flCity»New 

^ York, 1956. 

* • • • 

2 . Vonbun, F.O., “Tracking and Commwiications” pre«nted«| AlAA4ttl AntJUll * 

* • * • 

• • . 

meeting, paper 67*976, Anaheim, C«iifotnt»»Ociober 23-2T, 1967. 

3. Vonbun. F.O. and J. T. Mengel, “Tracking and Communications fo* Flanelary liarattid 

* • 

Missions”, Journal of Spacecraft and Rockets, Vol. 5, No. 7, pp 863<865. July 1968, ^ 

4. Martin, C.F., “Accuracy of Orbits Obtainable by Synchronous Satellites Tracking,** 
in Dynamics of Satellites ( 1 969), 1 20-1 29, edited by B. Morando (Berlin, Spring 
Verlag, 1970. 

5. Martin, C.F., “Evaluation of Satellites Tracking Satellites for Orbit Ot^ffminatioo and. 
Geopptential Recovery,” Final Technical Report on Contract NAS 5-t t7l6*Mod J, 
June 1970. 

6. •, CoQley, JJL.^and A. Marlow, “Orbital Errof Studies - Tracking frotn a SyndiroiMuft • 

Q • , 

Spacecraft,” GSFCX-55 1 -q^- 7, January 1969. 

9 V * 

7. Martin, C.F., T.V. Martin and D.E. Smith, “Satellite-to-Satellite Tracking fof EtftitMl* 

» » 

9 

ing Geopotential Coefficients.” The Use of Artificial Satellites forCe«dr»y, AGU 
Mono. 15, 139*144, 1972. 

• * • . . ^ 

' 8. - Schwar?, C.R., “Refinement of the Gravity Ticld by Safctlite*to*Salclliltf ‘ 

* < • * 

* 

tyacfctng.** T he-Dsc of Artificial Satellites for Geodesy, AGU Mono. IS, 133*138, 
.* 197 ?. \ 

. 9. • t‘Errof Earth PhyiicjiSai«^llilc Systems,” Final Rc|wl. Eart 

* 

. on NASA Cfaiil No. NCR 05-007*280. Ociol^ 1972. 

10. Hcf)'er»»n. P.J., “ATS4^/NIMBUS4F Ticking Data Relay Cxpcftmenl TciTiiUcat 
“ Sufbmary**.NASA<iSF'lC. JOctobet I968‘. 


1 1 . Vonbun, F.O., ct af., “Gravity Anon**ly - Apollo/Soyu/’^NASA^tSFC . 

X- 920-75-308, December, 1975. •' ‘ ' ' ’ . 

12. Teles, J., Ayres, C„ “Advanced Spacecraft Tlwking Techniqueix. .ing tlw Tricking 

and Data Relay Satellite System”, presented a| thdXAS/AIM AsttitdynamiciConfer- 

• ^ 

ence. Sept. 7-9, 1977, Grand Teton National Park, Wyoming. 

13. Schmid, P.E. and F.O. Vonbun, “The ATS-F/NIMBUS-F Tracking and Orbit Deter- 
mination Experiment”. 1974 IEEE INTERCON, New York City. 

* 

14. Schmid, P.E., F.O. Vonbun. and B.J. Trudell, “ATS-6 Satellite-to-Sa*tetCte Tracking 
and Data Relay Experiments”, pp. 1048-1058, IEEE Transactions on Aerospace and 
Electronic Systems, Vol. AES-1 1. No. 6, November 1975. 

15. Vonbun, F.O., "Geodetic Satellite Missions and GEOS-C Spacecraft”, Space Research 

XI- Akademie-Verlag, Berlin 1971. 

16. Vonbun, F.O. et al., Geodynamics Expetiment MA-1 28, Apollo-Soyuz test Project, 

. Preliminary Science Report. NASA TMX-58i73i, February 1976. 

17. Vonbun. F.O.,Tlte ATS-F/N1M8US-E Tweking Experiment”, pp. 112-120, 

’ ."Hota<ion of the Earth”, P. Mechipr & S. Yumi (eds), D. Reidel Publishing Co., 
tnirdrecht. Holland, 1972, 

1». Schmid, P.E.. "The Tra«,‘king. 4 ii«f Data Relay, Experiment ( f&DRE)”, pp. 207-218, 
The NHlBGS-6 User’s Guide NASA-<iSFC. February 1975. 

19. Marini, J.W., '‘Ooppler F4ctorsiii Satellil*-to-Satellite Tracking”, NAAA-GSFC 
X-93 2-74.93. April 1974. 

• • • 

2^.^ Schmid, P.E., R.B. Bent, S.K. Llfwelfyn, G. Nesterezuk, and S. Rangaswamy, NASA- 

« 

« 

'XlSFC lonosphwioj'orrectioni'to Satellite Tracking Data”, NASA-GSFC X-591-73- 
281, December 1973. 



r- 


21 . Marini, J.W., “Correction of Satellite Tracking Data for an Arbitrary Tropospheric 
Profile”, Radio-Science Volume 7, No. 2, pages 223-231, February 1972. 

22. Marini, J.W., “Tropospheric Range-Rate Tracking Data Correction”, NASA-CiSFC 
X-55 1-72-277, August 1972. 

23. Argentiero, P., and Loveless, F.. “Orbit Detemiination with the Tracking Data Relay 
Satellite System”, NASA-C.SFC X-932-76-185. February 1977 

24. Clark, George Q.. “Tracking and Data Relay Satellite Systems (TDRSS) User's Guide”, 
.STDN. No. 1 11.2. NASA-GSFCGreenbelt, Maryland, May 1975. 

25. Fang, B., and Gibbs, B.. “TDRSS Kra Orbit Detennii.ation System Kcview Study”. 
Planetary Sciences Department Report No. MTO 10-75. Wolf Research and l)evelo[v 
meiit Group. December 1975. 

26. .Argentiero. P., and Gar/.a-Robles. R . X'll OS-C Orbit Dctermin »tion With Satellue- 
ti>-Satellite Tracking”. Journal of the Astriinautical Sciences, Vo|. X.XllI, No. 3, pp 
241-256. July - September 1975. 

27. • Schmid. P.F.. Argentiero, P.. and Vonbun, F.O., “Satellite-to-Satellite System and 

Orbital Frroi i.stimates”, presented at the Precision Time and Time Interval Confet- 
ence, NASA-fiSFC. December 2-4, 1975. 

28. Schmid. P.1-. and J.J. Lynn. ‘Results of the 3 Novemlvr 1974 Appheatjons Technol- 
ogy Satellite-6 (ATS-6) Trilateration Test". NASA/GSI C X-932-75-104, April 1975. 

29. Lynn. J.J., P.L. Schmid and R.f. Anderson. "A New Method for Satellite Orbit 
Determination Using an Operational Worldwide I'ransponder Network". NAS.A-tJSFC 
X-59 1-74-2, January 1974. 

3C. Vonbun. F.O., “Satellite Trajectory Detetmination and Their I \pectej Firors OGO- 
IV, GEOS-I ” Dynamics of Satellites ( 1969), p. 89-1 03. Springer-Veilog. Berlin. 
Heidelberg. New Y<'rk. 1970. 

31 . Marsh, J., Goddard Space Flight Center, Private Communicatierns. 


152 


... _ , '^s 

N79-141 29 


ON-BOARD LANimARK NAVIGATION AND ATTITUDE REFERENCE 
PARALLEL PROCESSOR SYSTEM 

Lloyd E. Gilbert and Dhena T. Mahajan 
Martin Marietta Aerospace 
Denver, Colorado 


ABSTRACT 

An approach to Autonomous Navigation and Attitude Reference for 
Earth observing spacecraft is being developed. The technique Is to 
incorporate Landmark identification into the spacecraft on-board navi- 
gation and attitude control system. A fast landmark detection and 
registration system based upon a Sequential Similarity Detection 
Algorithm (SSDA) is being examined and laboratory experiments under- 
taken to determine if better than one pixel accuracy in registration 
can be achieved consistent with on-board processor timing and capacity 
constialnts. The SSDA is to be implemented using a multi-microprocessor 
system Including synchronization logic and chip library. The data is 
processed in parallel stages, effectively reducing the time to match 
the small known image within a larger image as seen by the on-board 
image system. Shared memory is incorporated in the system to help 
communicate intermediate results among microprocessors. The functions 
include finding mean values and suipiation of absolute differences over 
the image search rea Ine hardware is nlanned to be a low power, 
compact unit suitaule to on-board application with the flexibility to 
provide for different parameters depending upon the environment. 

INTRODUCTION TO LANIWARK TECHNIQUE 

The concept '<f using Landmarks to register images is common in the 
field of image processing. ^ Landmarks, also known as Ground Control 
Points (GCP), Registration Control Points (RCP) or anchor points are 
small Images with known geophysical coordinates. The known Landmark is 
found in a larger scene and thus the larger scene (at least the local 
area in the scene) is registered. The technique involves finding the 
best fit of a ’’chip" in a “window." A chip is a small image (size 
varies from 0x8 pixels to 32x32 pixels) of known latitude and longitude. 
The window is the large area to be searched, its size appropriate to 
the amount of uncertainty in where the chip will match. Figure 1 
shows an example of a chip/window pair. By finding the location of the 
chip in the window, the whole image can be registered. Many examples 
are found in the literature, for example Cloud Tracking from ATS 
pictures.^ 
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Known Latitude 
and Longitude 



Figure 1. Illustration of a Chip/window Pair 

The chip contains the knovm Landmark. The window is the area to be 
searched. 


The chip is compared to a possible location on the window by doing 
a one pixel at a time comparison over a placement of the chip on the 
window. The chip is then moved and the comparison r^'peated. A best fit 
is chosen. This gives a best whole-nixel match. Ta ic are two 
predominant styles of comparing chips to windows. The classical 
correlation coefficient Involving square roots of sums and products 
requires that the calculations be carried out over every pixel of a 
chip/window placement before a numerical answer is derived. Sequential 
Similarity Detection Algorithms^ involve sums of absolute differences 
between chip and window pixels and may be terminated before comparing 
every pixel of a chip/window placement.^ Using an SSDA approach with 
a dccroasinr threshold to allow onlv partial proccssii\c of most chip/ 
window placcn Ms^ a match can bt- quickly found. 
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After the best irtiole pixel match is found there arc several 
techniques to obtain best subpixel registration. These involve image 
enhancement^ or image resampling^; in any case interpolation between 
whole pixel placements is required. (It should be noted that there are 
several techniques to prepare the ir;:3ges for subpixcl evaluation. 
Although the raw data can be cuntpareJ, some techniques involve edge 
detection or contour following.) 

Implications on Spacecraft Navigation and Attitude Reference 

The idea of using Landmark data to determine spacecraft attitude 
and ephemeris information and incidentally register the images taken 
by the spacecraft is workable. Currently Goddard Space Flight Center 
uses the Landmark technique exclusively in their NAVPAK system. Tnis 
system completely registers images from Synchronous Meteorological 
Satellite (SMS). The NAVPAK output provides for the updating of orbit/ 
attitude state parameters exclusive of ranging or any other data. 7 
The Landmark technique has been shown to be better than traditional 
satellite tracking methods. 

• 

Various studies have been made which show not only the feasibi- 
lity of doing spacecraft attitude and orbit determination but have 
shown that the knowledge gained can be of higher accuracy than that 
derived from control system sensors. A 1971 work for SAMSO^ studied 
the mathematical techniques necessary to determine the attitude of a 
spinning geosynchronous satellite. Using as few as four Landmarks and 
a Kalman filter, the state vector of the spacecraft could be adequately 
described. The error analysis showed that high accuracy could be 
achieved in a few scans containing between two and four carefully 
selected Landmarks. 

A study of potential attitude and orbit determination and image 
registration techniques for the Earth Observatory Satellite^ examined 
various mixes of traditional and Landmark methods. Their study 
concluded that the attitude control system design must be considered 
as a part of the image positioning problem. The opposite is also true, 
the imaging system can be considered as part of the guidance and 
navigation system. The combination of onboard sensors (such as gyros) 
and Landmark identification can be illustrated to be a good control 
mechanism as well as enhancing the image processing procedure. There 
are various papers given at this symposium which address the problem 
of autonomous, on-board spacecraft calculations of attitude and orbit 
information. 

It is proposed at this time that the Landmark technique could be 
applied to a processor on-board a satellite to provide autonomous 
attitude and ephemeris update. 
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THE TECHNIQUE TO USE A SSDA TO 
ACHIEVE MEASUREMENTS FOR A GUIDANCE 
AND NAVIGATION SYSTEM 

Figurt; 2 Illustrates the general 
sequence of events. The Spacecraft (S/C) 
Is orbiting the earth with an imaging 
system which has a ground track such as 
shown in Figure 3. The image processing 
system contains numerous known Landmarks 
in its memory such as the chip shown in 
Figure 1. The estimated location of the 
chip in the field of view of the S/C can 
be calculated from the current S/C 
attitude and orbital knowledge. The 
location where the chip is actually 
found in the input data vs the 
estimated location generates error 
values wl>ich can be used to update 
attitude and ephemerus via techniques 
such as Kalman filters. 

The generation of this 
measurement on-board the spacecraft 
is the topic of this paper. The two 
algorithms of interest are the SSDA 
and the resampling approach. 



WAIT FOR S/C 
TO APPROACH A 
KNOWN LANDMARK 


>i 


OBTAIN IMAGE 
DATA; RESAMPLE 
DATA TO REMOVE 
CAMERA AND 
EARTH ROTATION 
SKEWS 

... 

USE SSDA TO 
LOCATE BEST 
ONE-TO-ONE 
MATCH OF 
INPUT TO 
KNOWN LANDMARK 


RESAMPLE (AT 
BEST MATCH 
LOCATION) TO 
LOCATE BEST 
0.1 PIXEL 
LOCATION 




- Knov^n Uiianyriis 


REPO'IT DELTA 
FROM T mated 

POS -N TO 
GUIDANa .AND 
NAVIGATION 
SYSTEM 


Figure 2. 


General Sequence 
of Events to Ob- 
tain a Landmark 
Data Point 


Figure 3. Image Sensor Ground Track 
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The SSDA Is traditional with the addition of a pseudo normaliza- 
tion to account for shifts In the mean value of the area of Interest 
due to camera gain changes or lighting conditions on the Earth's 
surface. The SSDA equation is: 


SSDA Value 


n n 

1 S Abs. Value ( (CPij -Chip Mean) - (WPij -Window Mean)) 
1*1 j“l 


where I 


implies sum over every pixel of the n x n Chip and the n x n 
area it is covering in the m x m Window 



is a pixel from the known chip 


WPij 


is the corresponding pixel in the unknown window at 
this placement 


Chip Mean * 2 £ Chip Pixels -r n 

i-1 j«l 


n n 

Window Mean “2 £ Window Pixels n^ for the current 

i»l j*l placement of the chip 


The chip is placed at a trial position of the window, the mean of 
the window under that position is taken, and the SSDA value for the 
sum of the absolute differences between chip and window on a pixei-by- 
pixel basis is determined. The chip is then moved to the next trial 
position and the process repeated. The best fit is the 
location where the SSDA value is a minimum. (A perfect match would 
result in an SSDA value of zero.) Note that the SSDA summation can 
be terminated when the summation exceeds any previous summation. 

After the best one-to-one match location is determined, the chip is 
resampled at 0.1 pixel intervals along-line and along-element axes. (Figure 4a.) 
The minimum SSDA value along each axis is the starting point for 
off-axis calculations. The nine locations surrounding the Intersection 
of the on-axis minimums are calculated as shown in Figure 4b. After 
this, values are generated to determine the minimum subpixel SSDA 
location. Figure 4b shows an example of two additional sets of 
measurements being required to surround the minimum value (*)• 
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ResanpK antf odtain SSOA aion^ eiements and aion9 roMs independently. 


Figure 4a 
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Resanole and oMair SSDA at 9 surrounding intersections of 
on ‘axis n^inimums. 

"Chase** minimum until smallest SSDA is located. 

Figure 4b 
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Several resample techniques are well known, among these are nearest 
neighbor (NN)» bilinear and cubic convolution (CC)> The nearest 
neighbor is useless in this application in that it simply duplicates 
the nearest pixel value and would not produce any change when attempting 
to generate a new sub-image. Bilinear is a straight line interpolation 
between adjacent values and has been found to be the best Interpolator 
for discontinuous data where higher order techniques tend to produce 
unjustifiably "undulating" values. For continuous data, that is, 
image data where adjacent pixels overlap or nearly overlap, higher 
order interpolation techniques such as the CC have been found to be 
highly successful interpolation schemes. The general form of a CC is: 


new old old old old 

where : 


‘^l - 

Cl - C 2 

(1 + dP) + 

C 3 (1 + dP )2 - (1 + dP)3 

K 2 - 


(dP)2 + 

(dP)3 


^5 ■ ^6 

(1 - dP )2 

+ (1 - dP)3 

•'a- 

Cl - C 2 

(2 - dP) + 

C 3 (2 - dP)2 - (2 - dP)3 


where dP * Subpixel displacement. 


Another technique has been found to be highly successful in chip 
placement. This technique will be called the "bilinear exaggerator" 

(BiEx). Before a new pixel can be generated, the slopes surrounding 
the area where the new pixel is to be generated are examined. If a 
trend is apparent which indicates the new value is not on the slope 
between the current pixel and its nearest neighbor then the modified 
slope of the preceding segment is extrapolated to generate the new 
pixel. (See Figure 5) An example of this case is: if the preceding 

pixel and the current pixel indicate a slope toward zero value and the next 
two subsequent pixels indicate a slope away from zero, then a local 
valley is indicated and the new pixel is generated based upon the slope 
from the preceding pixel. On the other hand, if there is a trend 
defined which indicates a continuous change, then a standard bilinear 
approach is used, i.e., the new pixel is generated on the slope between 
the current pixel and its neighbor. It should be noted that a new 
pixel will never be generated that is more than half a pixel distance 
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away. If the reaaople is required for greater than 0.5, the resample 
technique is run in the opposite direction. (Thus a move right of 
0.7 is accomplished by calculating a move left of 0.3.) This assures 
that the exaggerator will not produce unrealistic values. For a two 
dimensional resample the same technique is used with one exception. 

To generate a new pixel there exists a term which contains DL times DE, 
where OL is the move along a line and DE is the move along elements. 
This term is always treated as a bilinear case. (To continue to use 
the slopes in either the along-line or along-element direction would 
introduce a bias which has no physical justification.) 
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Figure 5. Resample Technique 
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t As chips are chosen which intent i%nall|Aontain 

Itfife (coast linest cross-roads, etc.) BlEx pvoves Highly^ 
succossl^l as it tends to exaggo^ito tho very discon^nuitdes 
chip was chosen for and thus generator, a "similarity v4th teiy 

sharp edges. Additionally, due to tne relative spe^s in ^igica^ 
processors, the selection of the proper slope for resampling %Ae«^ 
less time than the additional mult IplicaCion^tn the CC techn#<^. 
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EXPERIMENTAL RESULTS 
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The use of the SSDA to achlevoKubplxel accuracies has been 
examlt^jjLin a laboratory environment. A PDF computer fr^th a 

video ^Kplay is utilized. Programs have been written in Fortran 
with no attempt at optimization with respect to size or speech «f 
execution with one exception. The exception is that a moderadTe 
attempt has been made to stay with fixed point arlthi^tie where ^ve^t, 
possible. (The PDF 11/45 does fixed point arithmetic in 2 to 4 .44 s. 

For floating point an add takes apprc^mately 7 ms and • mult(j}ly 
approximately lO^t-s.) These programs accept manetlc tapes which 
contain image data from either Landsat or SMS spacecraft.?' The 
programs require extensive memory because oM the statistic Ic^eping - 
and reporting done in the laboratory environment. The actual J.and* 
marl: registration programs have been kept isolated however^' and requlr# 
apprcu||^tely 3000 words (g bits) of memory. Ihe whole pixel 
scarc^^equires approximately 750 words and the subpixed search v 

requires an additional 2300 words^^The al|^llty to use various re- 
sample techniques as well as a clais4ical correlation or SSIJA Whole 
pixel search technique is available. 

(o 

The program allows an operator to select fte size of a chip 
and the size of the window. Ihe following tables give some not 
typical timing requirements. 




■ 'a) 

( 


3 


3 


Table I, 


window Size 


|iole Pixel Search Timing Requirements 
Chip Size 


.is I, 


Time Required ^etonds) 


40 X 40 Pixels 

30 X 

30 

Pixels 

8.8 

40 X 4(jglxels 

26 X 

26 

Pixels 

12.3 

40 X 40 Pixels 

22 X 

22 

Pixels 

14.3 

40 x 40 Pixels 

18 X 

18 

Pixels 

14.3 

40 x 40 Pixels 

14 X 

14 

Pixels 

12.3 
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Table I. Whole Pixel Search Timing Requirements (Oont.) 

Window Size Chip Size Time Required (Seconds) 

40 X 40 Pixels 10 X 10 Pixels 8.8 

32 X 32 Pixels x 16 Pixels 3.6 

76 X 76 Pixels 16 x 16 Pixels 69.5 

Table II. Subpixel Resample and Search Timing Requirements 
Ctilp Size Time Required (Seconds) 


30 

X 

30 

Pixels 

19.0 

to 

27.0 

26 

X 

26 

Pixels 

17.0 

to 

18.0 

22 

X 

22 

Pixels 

10.0 

to 

13.0 

18 

X 

18 

Pixels 

8.5 

to 

10.0 

14 

X 

14 

Pixels 

4.5 

to 

6.P 

10 

X 

10 

Pixels 

2.5 

to 

3.0 


The times shown in Tables I and II reflect the experimental 
mode of operation where each SSDA Is summed to completion so that 
analysis may be done on the approach to the mlnlmums. The programs 
may also be run with the SSDA terminated whenever a previous minimum 
Is exceeded. Experiments show that a 6 to 1 Improvement In execution 
line is achieved In this mode. (In this mode a seed value is cal- 
culated based upon the mean value of the chip.) Additionally, 
analysis of the execution of the programs indicates that a 3 to 1 
Improvement in run time can be achieved If the resample algorithm 
Is linked to the SSDA such that only enough pixels are resampled to 
provide Input to a terminating SSDA. 

The results of numerous runs shows that the SSDA will whenever 
possible find a correct match when the subpixel algorithm Is Included. 
That Is to say, the whole pixel SSDA results occasionally point to 
an adjacent whole pixel location; however, the subpixel algorithm 
will then work Its way to the proper location. Thus a subpixel offset 
of 0.7 pixel Indicates a better match closer to the adjacent pixel 
than to the minimum valued one. 


* Exp*»|mentr Vtt1» «nrfe*ognizablc «hlps <h«avy cloud coyer OA 
wlndov tape where the ehtp wao ehostn frote a aloud- free tape) <.atio«r 

lJi.it llio SSD.I txhiblts a charnrtt ::istic waivJcrinn xcJicn no match i* 
availahffc» Calculations 4n the j r '/.ram oau.io a "N’o-lUitch'* indicator 
to be*scC urfJe# these conditions. I'.sporimcnt# witl> cliips Which were 
'jpuappscly chosen to he of quastiotnh^e quality ^P^rtlal cloud rover* 
tfiie ^imenSiona) chartctcristlcs) indicate that the SSDA is senatCiv* 
to slight mis/registirattons and has characteristics which atiow a 
"tjualtty** value to be placed upon any given fnatch. 

• 

figure 6 shows typical SSDA responses, to an excellent rtiip 
(clear^ two dimensional characteristics) a poor chip (partial stood 
cover o« chip or window or one dimensional characteristics) and art 
inadequate chip (heavy cloud cover ot a nondescript chip •clvdt|on)« 

ExpertmenCal. Conclusions 

While we ate sglll in the process of developing s«ch expcriMmfat 
knowledge as the correct size of the chip and the optimum retanvllng 
techniques* some Cenative conclusions caa be reached* The three 
most Important ate: 

The gSDA wilt generate a minimum Walue at the correct location 
to + 0.1 pixel with continuous Input data*. (However* note that the 
effect of all possible noise sources hps not yet been Included In the 
laboratory experiments.) 

The SSDA will generate a unique profile that indicates the • 
quality of the match. 

The search pattern for the best subpixel toeaCion works ttt alt 
cases; there is no need to generate every possible subplxet value* 

IMPLEMENTATIOM 

Overview 

> 

TTie Sequential Slmllail^lty Detection Algorithm (SSDA) descrlbc4 
In the previous section lends itself for parsllel processing# A 
mult-microprocessor system (MMS) 4s props»r4 to implement tfkc SSDA« 

Simple processing mlements vith only moderate processing power #r« 
proposed Aestus# of the elementary Aature of the eonputatione involved^ 
The MMS is well euited for space-borne appUcattonsi 

Vslng a serial conventional computer to Implement the SSDA poses 
three main problems; weighty power and space* Although most Conventional 
computers have a large processing power ^ they arc not cuCtable tof Space* 
borne applications due to thsit Mgh volume^ weight and power eon* 
sumption# It is not practical, from the point of View of real ttne 








response, co use single processor, either. Iherefors the suthists V 
have proposed a light 'flight, low power mult-mlcroprocessor syskein« 

It l3 vstinatcJ that tlu M'!S will he able tv> process the data rfNfelvod 
irvxr the ltaa.;e System an I provide response to Gtiidance and Contl‘<»l * 
System almost in real tine. 

Hie main point of this section is to show that with currently 
available microprocessors and RAM memory devices a system can be 
^ economically developed for Autonomous Landmark Navagatlon and 
Attitude Reference task on-board rather than using ground support 
computer. An estimate is made of the time and cost required Co solve 
this problem on the proposed system and compared with the simulation 
values obtained for a DEC PDF 11/45. 

Ihe following sections describe the MMS salient hardware and 
software features. It is a system tailored to a specific application, 
i.e. Implementation of SSDA, and no generality is Intended. Ihe 
overall system diagram Is shown below, in Figure 7. 



Figure 7. Overall System Diagram 


The scope of this section is limited to describing the processlt\K' 
function of the MMS. Other necessary features for a space-home 
application such gs fault-Culerance, ruggedized design and down- 
loading communications . Interface are reserved for future study. 

Architecture 


The following set-up is proposed for on-board implementation of the 
SSDA, using the MMS, Guidance and Control System and Image System. 





Figure System Block Diagram 

Ihe MMS comprises of Processing Unit (FU), Synchronization Logic 
(SL) end Chip Library (CL). Figure 8 shows System Block Diagram with 
Interconnection •f the MNS elements to the Image System and Guidance 
& Control Systei 9 «> Ihis paper assumes that the following inputs are 
available to tils MMS. 4 

# 

' 9 ' 

X'aage 43 l*tem w Die Image system views a strip of 240 Km width on 
Q the «arth surface* Die picture is quantized into 

O gfay scale values. Die scanning mechanism scans 

^ along the width of the strip and outputs 6 pixel 
values In parallel. One scan contains (6 x 3000) 
o pixels and takes about 70 nilli-seconds. Die pixel 

0 0 values are eight bit quantities and vary from X'OO' 

0 ^ through X'FF*. Here X denotes hexadecimal values, 

oo the Image system also sends a signal to Indicate 

^ 0 ° o start of a scan, and a signal to Indicate transmlfslon 

0 of a pixel. 


£? 


o 


& C System 
o 0 
o 0 


o 


& 


Die Guidance, Control of Navigation System outputs 
coordinates Indicating when in a scan the MMS should 
start collecting Window pixel data received from 
<the Image System. Die G & C System also outputs 
Chip Select information so that the MIdS can pull 
the chip data from library and store it in memory 
modoles Ml through M 6 . 

Input to the G & C System is a set of coordinates 




(x« y) Indicating where in window the best fit 
occured for a given chip* These coordinates are with 
respect t.> the window position outputed by the G & C 
System. The G & C System, then, may take any 
necessary navigation and attitude reference corrective 
action. 


Synchronization 

A synchronization logic is a unit of the MMS* that interfaces 
to the G & C System and the Image System. Its function is fully 
implemented in hardware, n>slnly because of the simplicity and 
invariance of the function. 


It is assumed that the G & C System will send the followik.g data 
to locate a window 


^I't Number of scans to elapse from now 
( 2 ) Pixel number in the scan of interest where 
the window starts. 


From 

Image System 


I 

V 





Data from 





& C System 


J 

f 



a Clock 1 

# of Scans 

H 

Pixel # in Scan 


Pixel Clock 


Figure 9. G and 0 System 

The G & C System will set up two registers in the synchronization 
logic as shown in Figure 9. These registers are simple count down 
registers using appropriate clocks from the Image Sys*"em as shown in 
the above figure. When the "Nfumber of Scans" register reduces to 
zero, the "Pixel number in scar" register starts counting down. When 
the later reduces to zero a "start" control f Ignal is sent to the PU. 


Chip Library 

The (Slip Library unit of the NMS holds pixel data for chips. 

The number of chips should be sufficient for Landmark Navigation and 
Attitude Refzrerce task all over the surface of the earth. The authors 


167 





propose to implement the CSiip Library using bubble memory technology, 
because It is a more reliable and compact auxiliary memory unit as 
compared to the error prone magnetic tape unit. Table III below shows 
perforraince of a TI magnetic-bubble memory system. 

Table III 

92 Kilobytes 
0.69 lb 
38 in3 
44 kb/s 
11.5 Watts 

10-9nd per bit n = years residence 

time 

d * operating duty 
cycle 

One 92 Kilobyte memory system is sufficient for 92 30 x 30 chip 
or for 276 18 x 18 chips. Further study is required to determine size 
and number of chips necessary for the Navigation and Attitude Reference 
problem. 

Processing Unit 

The Processing Unit is the most iinportant unit of the MMS. 

Figure 10 shows Organization of the Processing Unit. There are two 
types of r.ain components: Processing Elements (PE) and Memory 

Modules (M)* The PEs are arranged to process parallel data received 
from the Image System and from the Chip Library. They are also 
arranged to form a pipeline to process received data. The Memory 
M.-:dules are used to hold and transfer data from one processor to 
another. 

Processing Element: A Processing Element consists of a 8-bit, fixed- 

tnst.uction-set microprocessor. We have chosen 8-bit wide microprocessor 
beca the pixels are 8-bit wide. However, the microprocessor must 
have instructions to manipulate 16-bit data because the SSDA values 
are 16-bit wide. The third requirement is to have at least 2 incex 
registers to effectively manipulate window and chip data, each stored 
in a matrix form. There are number of microprocessors currently 
available and the authors think that Zilog Z-80 microprocessor 
suits well for the Processing Element. We have chosen a fixed- 
instruction-set microprocessor as against a mlcroprograninable 
processor because the MMS is to perform a well defined dedicated 
function rather than to be used for general purpose processing tasks. 

The choice also reduces space requirement. Each of the Processing 


Capacity 

Weight 

Volume 

Transfer rate 
Power dissipation 
Ha’-d error rate 
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end: 



Figure 10. Processing Unit 













Elements PEI : PE6 has two input ports. One input port is to get 
chip pixel data from the Chip Library unit. Other input port is to 
get window pixel data from Image System. Ihe Processing Elements 
PE7 : PE9 do nor have any I/O port. The Processing Element ElO has 
one output port to send best fit coordinates to the G & C System. 

Each of the microprocessors has an address space of 64 K bytes. This 
space is divided into three categories: Read only Memory (ROM)^ 

local memory and shared memory. The ROM holds program to handle 
power up, down loading etc. The local memory holds part of a 
program peculiar to the microprocessor to implement the SSDA algorithm. 
The local memory also holds. local data. 

Shared Memory 

The shared memory holds data that is to be transferred from one 
processor to another. The shared memory modules are labled Ml 
through M7 in Figure 10. The details are shown in Figure 11. 



Figure 11. Dual Port Shared Memory Module 


A shared memory module has two (for Ml : M6) or four (M7) ports 
connected to the microprocessors. We describe operation of a 2- 
port memory which is similar to 4-port memory module. TWo ports 
are controlled by clocks derived from a single clock as shown 
below. 


Common 

Clock 



1 Clock 

2 Clock 
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Memory access requests from microprocessors are sampled at the, say, 
rising edge of the clock. Once a request Is detected from one micro- 
processor, requests from other microprocessor are blocked by sending 
WAIT signal to that microprocessor. Once the request from the first 
microprocessor is completed, secone microprocessor Is allowed to 
access the memory. This technique assures mutual exclusiveness of 
memory accesses in shared memory module from two microprocessors. 

An example is shown below to illustrate the point. 

Let ABRF ■ address of flag in shared memory 
ADRD ■ address of data In shared memory 

Assume that the flag is reset. Data is to be passed from PEI to PE7. 
Assuming Zilog 2-80 microprocessor instruction set, the following 
codes accomplish the cask. 


PEI 


PE7 


LOOPl LD 

(ADKF) 

LOOP? LD 

(ADRF) 

ADD 

A 

ADD 

A 

JP 

NZ, LOOPl 

JP 

Z, LOOP 

LD 

(ADRD), HL 

LD 

HL, (ADRD) 

LD 

A, 1 

LD 

A, 0 

LD 

(ADRF), A 

LD 

(ADRF), A 


In the example, the loop at LOOPl assures PEI that previous 
data is processed by PE7. Ihen PEI stores data to be transferred 
at (ADRD) and sets a flag at (ADRF). The loop at LOOP? assures PE7 
that the data to be transfered is available. Then PE7 gets that data 
and resets a flag to Indicate to PEI that next data can be transferred. 

Parallel Processing 

We consider some programming aspects in this section. We describe a 
few s'eps of the SSDA algorithm to uhow the parallel nature of the 
processing involved. Two processing steps. Input and SSDA computation, 
are shown in Figures 12 and Figure 13 respectively. . Horizontal lines 
indicate process, whereas blanks indicate idle time. Figure 12 shows 
that the process PEI : PE 6 starts after receiving a signal from the 
Synchronization Logic. AIT 6 microprocessors accept data from the 
Image System in parallel and store it in their local memories. Each 
microprocessor is programmed to accept the number of pixels equal to the 
length of the window along the scan direction. Then each processor 
computes partial sums from the data Just received. If we have a m x m 
window and nxn chip, then (m-n+1) partial sums have to be computed 
per scan. Then the microprocessors wait till a "continue" signal 
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Flmre 13. Parallel Proce<;sing : SSDA Computation 



Is received from the Synchronization Logic* Hie microprocessors 
repeat same operation until scan is completed* 

IXiriug the Input step, microprocessors P7 : PIO arc idle. 

However, they arc active during tiio SSil\ computation a:id subsequent 
steps. Figure 13 she ^ SSDA computation step split into 7 separate 
processes, A through G. We assume that the top of window corresponds 
to top of the first scan. To compute a SSDA value, the first step 
is to find a Window Mean, l.e. to add (nxn) window pixels at current 
coordinates and divide by n^. Adi operation reduces to adding n 
partial sums that were computed during the previous step^ For each 
microprocessor the operation further reduces to adding ( 7 ) partial 
sums. In the next process, B, microprocessors E7 : E9 add numbers 
supplied by El : E 6 and pass on 3 numbers into M7. In the third 
process, C, microprocessor ElO adds 3 numbers as they become available 
and divides the sum by n^. Since it is an unsigned division and over- 
flow, underflow conditions are ruled out, the division algorithm is 
simple, and takes less than 200yu s for 4 MHz Z-80 microprocessor. The 
Window Mean obtained in ElO is passed on to El : £6 by £7 : £9 in 
process D. 

Once the current Window Mean value is available in El : E 6 , 
a plpellne-like process is started. Hie microprocessors El: E 6 
compute the ABS (W-OfK) value, E7 : E9 add 2 numbers and pass on sum 
to ElO which accumulates sums to form a SSDA value for the current 
coordinates. 

PERFORMANCE PREDICTIONS 


Although it is almost impossible to analyse complete performance 
of the ^D 1 S here, we give below results of some timing calculations 
pertaining to two SSDA algorithm steps described in the Parallel 
Processing section. For our calculations, we have assumed a 4 MHz, 
Zilog Z-80 microprocessor to be used in a Processing Element. We 
have also assumed window and chip sizes to get some figures from the 
formula. 


Window 


* ra X m • 

Size 90 X 90 pixels 


■nxn® 

Chip size 30 x 30 pixels 


Input 

Image System: Sends 6 pixels (bytes) in parallel 

Total bytes sent • (6 x 3000) in 70 ms. 
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'C 


# of scans of Interest ■ ro ■ 90 °* 15 

6 6 

Input time *= m x 70 * 15 x 70 = 1050 ms. 

6 

MMS: 6 microprocessors accept data in parallel at the Image 

(El : E6) System rate of about 20 vls per pixel, although they 

can accept at a faster rate of lOi^ s per pixel. Input 
time per scan " m x 20 * 90 x 20 ti s • 1.8 ms 
Partial sums per scan ■ m-n+1 • 61 
Time to find partial sums ■ (m-p+1) x 25^<s 1*5 ms 
Idle time per scan ■ 70 - 1«8 - 1*5 • 66»7 ms 


Window Mean 

El : E6 Each processor adds (^) * 5 partial sums 
Time required •* 150.. s, 

E7 : E9 Add 2 numbers & pass on sum * 25 a« s 
ElO Add 3 numbers & divide * 250.. s 


E7 : E9 Pass on window mean to Ml : M6 • 15 s 


SSDA Computation 


El : E6 


E7 : E9 


To find one ABS (W-C+U) value “ 50/is 
Total time = n x n x 50 * 30 x 30 x 50 i.s * 7*5 ms 
6 6 

Add 2 ABS values and store it in M7 
No extra time due to pipeline effect. 


ElO Accumulate sums of ABS values, to form SSDA value, 

for current coordlitates. No extra time due to 
pipeline effect. 


Total time to find coarse location of best fit is equal to 7*5 x 
(m-n+1)^ « 7*5 x 3600 « 27*0 seconds which compares favorably with 70 
seconds obtained on PDF 11/45 for a Fortran program written for 76 x 76 
window and 16 x 16 chip. 


We list below the main components needed to build hardware for 
proposed multi-microprocessor system. The cost given Is approximate 
and does not include hardware/software development efforts. 

For a 90 X 90 pixel window and 30 x 30 pixel chip, the MMS needs 
following memory capacity. 
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Local Memory : 

PEt through PE6 

m 

4 

Kbytes 

each 


PE7 through fE9 

• 

1 

kbytes 

each 


PEIO 

m 

8 

kbytes 

each 

Shared Memory : 

Ml through M6 

m 

1 

Kbytes 

each 


M7 

a 

4 

Kbytes 

each 


Total Memory size * 45 Kbytes 

Using a 1 X 1 k military temperature raage RAM chips at about 
$8 each (including addressing logic) gives us an estimate of $9600. 

Cost of 10 microprocessors and support chips at $300 each gives 
us an estimate of $3000. 

Early price of one bubble moaiory device (TI 92 k blts)^^ was 
quoted at $200 and is expected to drop in 1978. 

The power consumption is expected to be of the order of 70 -75 
watts. The approximate break up is 15 - 20 watts for microprocessors, 
20 watts for RAM, 15 watts for bubble memory, and 20 watts for other 
logic chips. 

The weight of the MMS is estimated to be of the order of 20 -25 
lbs. including bubble memory and chasis. 




a 
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AUTONOMOUS SATELLITE ORBITAL 
NAVIGATION AND ATTITUDE DETERMINATION 

Shing Peter Kau 

Avionics Division 
Honeywell Inc. 


ABSTRACT 


A known linear landmark navigation system is described in this 
paper. It involves the use of an electro-optical sensor to 
provide sightings to linear earth features such as highways '' 

and coast lines. The sensor concept and the navigation system 
mechanization are described. Performance analysis results 
show that landmark sightings provide accurate navigation up- 
date and that this accuracy can be preserved using radar alti-* 
meter measurements. , ’ 

Description-on a stellar inertial attitude determination system 
is also presented. Attitude reference performance consistent 
with the requirement of the navigation system is shown to be . 

achievable by this method. • • 

■ * • * ’’ 

1 . 0 INTRODUCTION ^ * . • * 

This paper describes methods of autonomous satellite navigation;.. ‘ ' 
and attitude determination using on-board sensing and processing 
capabilities. Sensor concepts, system mechanization approaches, 
and projections of navigation and attitude reference performances 
are presented. The paper is divided into two parts for separated 
discussions on the navigation and attitude determinatior ■ “blems. 

2 . 0 AUTONOMOUS NAVIGATION VIA LANDMARK SIGHTINfiS 


Satellite navigation using knewn or unknown earth landi . . htv 

bean intensively explored in the past. References 1, 2 
Navigation information is derived from tracking known or cr" 
landmarks. This implies that the landmark sensor will be ' i 

with large o« movable Field of Vi -w (POV) and that the ian;' 

«ill be a point target or small area with distinct signatu.. . jving 
% veil defined centroid. For low altitude orbital application^', 
known landmark navigation approaches typically are less sensitive . 

to pointing errors than the unknown landmark approaches. *equire- * 

ments on attitude reference and landmark sensor accuracy are there- 
fore less stringent. The unknown landmark approaches, on the othef 
hand, are attractive in that the task of landr..ark idenlif ication 
can be eliminated, thus relieving the requirements on storage of 
landmark signillares. 0 o 
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Tlt« autonomous navigation mathod considered here uses such known 
linear earth features including highways and coast lines as * 

candidate landmarks. Due to their long physical dimension, strap- 
down sensor with relatively small POV can be m jhanized for de- 
tection of landmark crossings. Being a known landmark approach 
in nature, the proposed method has the advantage of low sensitA* 
vity to pointing errors. Thus, the design requirements for tlM 
landmark sensor ate quite relaxed. Also, the simple signatur# 
of linear features does not require extensive storage for known 
landmark catalog. 

2 . 1 SENSOR CONCEPT 


The landmark sensing down sensor is an electro-optical devies* 

It consists of a telescope that earth surface features 

onto two linear silicon detes:tor arra/s which arc separated by 
3 degrees and oriented 45 degrees with . ecpect to v.he direction 
of nominal image motion as depicted in Figure 1. Due to the cross 
array component of image nx^tion, the FOV of detectors over terrain 
scenery sidesteps from scan to scan. Hence, two dimensional dis- 
crete images of terrain scenery can be created form successive 
samplings of detector cell readouts- T^’.^se digital images are 
processed to derive landmark sighting information for system 
navigation update. 

The measurement provided by the down sensor is the LOS-vector to 
the centroid of the segment of a linear earth feature that falls 
into the sensor POV. This is obtained from processing the dis- 
crete image for detection of the presence of a linear feature 
and for extraction of the feature oriental- ion and the segment 
centroid location. Due to the deterministic signature of linear 
landmarks, deterministic image processing techniques such as 
thresholding and e<<ie enhancement are used. The data processing 
techniques and potential sensor accuracy will be dem >nstrated 
through the discussion of simulation data for a test case. 

The test case involves viewing a road of 50 ft width from 100 n.mi 
alti'ude under a 45 degree sun angle lighting condition, fhe 
linear feature is characterized by pavement with uniform reflect- 
ivity of 0.5. The background is represented by an exponentially 
correlated spatial rr<^-;ess with correlation distance of 500 ft, 
mean reflectivity of 0.25, and standard deviation of variation 
of reflectivity of 0.08. A poi'tion of the simulated original 
scene is shown in Figure 2a. with each letter representing the 
reflectivity, in st«ps of 0.1, cf the elementary area (:;20 x 20 ft^) . 
The down sensor detector cell width, scaled to 37 arc seconda, ha^ 
a ground projection of approximately 110 ft. With an array scan 
rage of 1000 Hz, the FOV of a detecto'' in consecutive scana has 
an overlap of 5/6 of a cell width. Th; s overlap of FOV, together 
with the multilevel cell readout, allows a limited degree of 
improvement of image resolution. A super-positioned imaqe created' 
from consecutive scans of the array cell readouts is shown in 
Figure 2b. The terrain radiometry, detector and electronics noises, 
and the 3-bit quantization oft cell readouts have been fullv siru- 
The result obtained from a simp* 3 thresholding operation 
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is rSown in Figure I'c. Finally, a directional jradient operation 
is . '.'formed on the resulting image shown in 2c to yield Figure 
2d with the linear feature significantly enhanced. The direct- 
ional gradient is defined as the dot product between a prespeci- 
fied unit vector S and the gradient vector of the image func- 
tion W(X,Y) . With W (X,Y) treated as mass den sities, the center 
of mass, the moment of inertias and the principal axes are com- 
puted. To facilitate the detection of the presence of linear 
feature, a shape factor is computed ap the ratio of moment of 
inertias about the principal axes. The center of mass provides 
the information of the location of the centroid of the linear 
feature segment. The orientation of the linear feature is pro- 
vided by the principal axes. A summary of the results obtained 
from the feature detection ard extraction processing on the three 
images generated from the intermediate steps of enhancenent is 
presented in Table 1. 

2 .2 AUTONOMOUS NAVIGAT1 Q> / SYSTEM MECtlANIZATION APPROACH 

The autonomous navigation system concept described in this paper 
consists of a down sensor for landmark sighting, an on-board 
computer and resident software, a long term stable clock, a radar 
altimeter, and an attitude reference subsystem. A functional 
block diagram description of the system concept is shown in Figure 
3. The nominal navigation solution is computed through the inte- 
gration of vehicle equatlo i of motion using modeled accelerations 
(drag, gravity). Due to errors in the initial conditions and 
uncertainties in the acceleration models, the error buildup of 
the nominal navigation solution requires periodic upviates using 
sensor measurements. These include the down sensor )*nown land- 
mark crossing detection and the altimeter altitude measurements. 
The optimal implementation of these measurements calls, naturally, 
for the application of Kalman filtering techniques. An im.portant 
operation involved in any approach of known landmark navigation 
is that oi landmark identification. The system mechanization 
approach will be outlined in the following in terms of the inter- 
pretation of dow~i sensor measurement geometry required in the 
Kalman filter formulation for implementing these measurements 
and the data processing flow of the landmark identification pro- 
cedure. 


2.2.1 Kalman Filter Form »' lation 

The measurement provided by down sense is the LOS-vector, denoted 
as ^ , to tha center of the segment of the linear feature that falls 
into the sensor FOV. Since the sensor FOV is small compared with 
the length of a linear lancr .rk, the exact point where the sensor 
LOS intercepts the linear feature is ambiguous. To circumvent 
chis ambiguity, a miss distance is computed from the LOS-vector 
measurement prior to Kalman filter navigation update procer.sing . 

Let the vehicle position vector be Sv *rid the intercept of L with 
the linear landmark be X (the target) as depicted in Figure 4. 

The targst position can be computed as: 



CT*J= (1) 

where the notauion convention is such that ( denotes the 
inertia] -frame coordinates of a given vector X . The slant 
range can be computed as: "" 

ss= , 2 , 

whPi i Ylg is earth radius » superscript T denotes matrix 
transpose and 

Let ^ be the urit vector normal to the plane containing the 
linear landmark* The miss distance between the projected down 
sensor target point and the linear landmark is computed as the 
dot product r 

In the case where perfect knowledge of vehicle position and 
LOS-vector were involved in the evaluation of equations (1) , 

(2) and (3)^ the resulting miss distance would be identically 
zero. The actual value of dot product reflects errors existed 
in the a priori knowledge of vehicle position , attitude, and 
down sensor LOS-vector measurements. A formal expression that 
relates miss distance to errors in navigation attitude 

( 458 ) and sensor errors is given as follows: 

d=H,(AlvJ+H^C«2®J +H3(5) 

Actual expressions for Hi , H 2 and H 3 can be obtained by differ- 
entiations of equations (1) , ( 2 ) and (3) . 

2.2.2 Landmark Identification 


The linear landmark involved in a down sensor measurement must 
be correctly identified to enable the extraction of useful navi- 
gation information from landmark sightings. A landmark catalog 
will be carried on-board to facilitate the identification pro- 
cedure. Each linear landmark is defined within the catalog in 
ter:.;s of the location of a reference point and the orientation 
of the feature with respect to local north. 

Upon a down sensor linear feature detection, candidate landmarks 
in the vicinity of the projected sensor FOV will be tested in 
two steps for identification. First, the orientation of candi- 
date landmarks will be compared against the measured feature 
orientation using a threshold established from expected attitude 
reference and feature orientation measurement errors. Miss 
distances to the candidate linear landmarks that survive the 
orientation screening are then c mputed. The miss distance, in 
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general, consists of two components. First, the actual miss 
from sensor target point to the candidate landmark. Second, 
the equivalent miss contributed by errors in the a priori 
navigation and attitude information and the sensor measurements. 
For the correct candidate, the first component is identically 
zero. 

The level of the second component can be predicted from the 
navigation and attitude covariance matrix evaluated as part of 
the Kalman filter computations. By careful selection of candi- 
date landmarks to avoid sightings from congested areas, the 
first component of the miss distance to incorrect candidate 
can be made considerably larger than the second component. This 
permits the discrimination between the correct and the incorrect 
candidates. A reasonableness test on the miss distance can be 
devised using a threshold computed from the statistics of a 
priori errors. 

The down sensor landmark measurement will be implemented for 
navigation if, and only if, a unique candidate is identified. 

A flow diagram summary of the landmark identification procedure 
is contained in Figure 5. 

2.3 NAVIGATION PERFORMANCE ANALYSIS 


The known linear landmark navigation performance results pre- 
sented in this section were obtained assuming a system that 
employed only one down sensor looking along the yaw axis of a 
local vertically stabilized vehicle. The down sensor errors are 
characterized by 0.15 mrad white noise and 0.15 mrad bias. The 
assumed attitude reference errors consisted of 0.15 mrad random 
noise and 0.15 mrad bias. A 9 state Kalman filter was assumed 
for navigation update using down sensor landmark sightings and 
radar altimeter measurement of vehicle vertical position. State 
variables considered in the filter formulation included 3 position, 
3 velocity and 3 residue errors for acceleration modelings. The 
ground track of the reference orbit used for performance analysis 
is plotted in Figure 6 indicating the down sensor landmark sight- 
ing schedule and the times altimeter is activated. The altimeter 
is activated only over ocean where the mean geoid height can be 
accurately modeled. Figure 7 contains the plots of the 3 axes 
RSS position errors from two covariance analysis runs. The curve 
labeled (a) is obtained implementing the landmark sightings 
scheduled for all two orbits. The curve labeled (b) is obtained 
implementing only the landmarks encountered in the first orbit. 
These results show that accurate navigation information can be 
derived from down sensor landmark measurements. Also, the navi- 
gation accuracy can be preserved over an extended period of 
orbital flights without additional landmark sightings. An impli- 
cation of this is that the frequency of landmark sightings required 
for high accuracy autonomous navigation can be rather small. Sys- 
tem errors assumed in generating these performance results are 
summarized in Table 2. 
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3-0 


AUTONOMOUS ATTITUDE DETERMINATION VIA STAR SIGHTINGS 


Satellite attitude determination can be accomplished by the inte- 
gration of the vehicle rates measured by a set of body fixed 
gyros. However, this attitude solution will diverge with time 
due to gyro drift errors. Periodic stellar updates can remove 
this attitude error build up and result in a high accuracy systam 
suitable for applications with extended mission duration. 

The stellar-inertial system considered here consists of a body 
mount star sensor, a set of three nominally orthogonal gyros, 
and an on-board computer. The star sensor consists of a teles- 
cope with a set of six detector slits placed on its image plane 
as depicted in Figure 8. A transit pulse is produced when the 
image of a star moves across a detecting slit. The basic star 
sensor measurement is the precise time when the star transit 
occurs. Kalman filtering technique is employed for optimal imple- 
mentatioi' of stellar attitude update. The descriptions of the 
update mechanization and the performance analysis results are 
presented in the following paragraphs. Th3se r suits are included 
to support the attitude reference system error budgeted in pre- 
ceding analysis of autonomous navigation performance. 

3.1 STELLAR UPDATE MECHANIZATION 


The autonomous attitude determination approach considered here 
involves the solution of vehicle nominal attitude and the update 
of this nominal attitude with star sensor transit time measure- 
ments. 


From the geometry shown in Figure 8, the condition for the star 
transit is described by the equation: 

d = N . S 0 


Where S is the LOS-vector to the transit producing star, N is 
a unit”"vector normal to the plane containing the detecting slit. 


To implement the transit measurement for attitude update, the 
dot products for all candidate stars and slits are evaluated 
using the nominal attitude 'Tax 


where 



Tr 


Tt ^ 


t h 

normal vector for i detector slit, 
LOS to candidate star 
measured transit time. 


(5) 
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The star and slit identification is accomplished by comparing 
all the dot products for a reasonableness test. Similar to 
previous discussions on landmark identification/ the dot product 
here consists ot two components. The first component is the 
angular miss distance between the candidate star and slit. This 
is identically zero for true transit producing star and the true 
detecting slit. The second component is due to errors in the 
nominal attitude solution and in the detected time of star tran- 
sit. 

The candidate dot products are compared with a threshold computed 
based upon the covariance matrix evaluated in the attitude Kalman 
filtering computations. The success of this identification pro- 
cedure lies in the fact that attitude errors are small when com- 
pared with angular spacing between detectable stars. 

Here/ again, the on-b»''ard star catalog can be tailored to avoid 
congested regions in tre celestial sphere. The star transit 
will be implemented for attitude update if, and only if, the 
star and the slit are uniauely identified. 

The dot product evaluated for the true transiting star and detect- 
ing slit provides a scalar measurement of the error in the nominal 
attitude solution. Through first order perturbation of equation 
(5), the value of dot (perturbed from the true value of zero) is 
related to attitude errors as; 

<i ( I'J ^ (0 i „ 

where 


CV) 


= roll, pitch, yaw attitude errors. 



O 


184 


'4 



Eauation (6) defines the :ittitude error observation provided by 
star sensor transit measu^ anents. A Kalman filter can be formu- 
lated on the basis of this equation to implement the sensor for 
attitude update. A block diagram showing the stellar attitude 
update processing is shown in Figure 9. Details on the develop- 
ment of the filtering equations can be found in Reference 4. 

3.2 ATTITUDE REFERENCE PERFORMANCE ANALYSIS 


The stellar inertial attitude reference performance results pre- 
sented in this section were obtained assuming a system that an- 
ployed only one star sensor with its LOS pointed outward 30® off 
the vehicle pitch plane. The star sensor errors are character- 
ized by 0.05 mrad white noise and 0.05 mrad bias srror. A 6 
state (3 attitude and 3 gyro bias) Kalman filter wai’ assumed for 
attitude update using the star sensor transit meassurements. 

Figure 10 contains the plots of the 3 axes RSS attitude errors 
from two covariance analysis runs. The curve labeled (a) is 
obtained assuming a local vertically stabilized vehicle. The 
curve labeled (b) is obtained assuming an attitude maneuver for 
star transit acquisition. The convergence characteristics is 
significantly improved. This maneuver, depicted in Figure 11, 
is designed to acquire complete observability on the six attitude 
and gyro bias states. The observability analysis leading to the 
selection of this maneuver can be found from Reference 4. Numer- 
ical values assumed in these analyses for various attitude ref- 
erence system error sources are summarized in Table 3. Notice 
that the attitude reference performances presented in Figure 10 
are consistent with that allocated in the navigation analysis. 

4 . 0 SUMMARY AND CONCLUSIONS 


Methods for autonomous satellite navigation using known linear 
landmark sightings and attitude determination using stellar- 
inertial sensor measurements have been presented in the above 
discussions. Performance analysis results obtained for the pro- 
posed autonomous navigation approach show that sightings to 
linear landmarks provide highly accurate navigation updates. 
Also, it is shown that the navigation accuracy can be preserved 
over extended periods of landmark free orbital flights using 
radar altimeter measurements. Frequency of landmark sightings 
necessary to satisfy given navigation performance goals can thus 
be relieved. Performance analysis results obtained for the 
stellar-inertial attitude reference system show that accuracy 
consistent with that required by the autonomous navigation. 




Figure 1. Dowti Sensor and Discrete Image of Terrain Scene 
Created from Detector Cell Readout Samples 
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Table !♦ Linear Landmark Measurement Performance 
for the Example Case 
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NO GRADIENT 

57 ARCGEC 

64 ARCGEC 

2.0 PEG 
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FIGURE 2 (B) 
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FIGURE 2 (C) 
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FIGURE 2 (0) 







Figure 3. Functional Block Diagram of 
Autonomous Navigation System 


187 
































Figure 4. Geometry of Linear Landmark 
Sighting 
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Figure 7. Navigation Error History Obtained 
from Covariance Analysis 
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Table 2. Nominal Values of Navigation 


ERROR SOURCES 

Error Sources 

STANDARD DEVIATION 

CORRELATION 

I. 

SENSORS 

DOWN SENSOR: RANDOM 

0.15 M RAD 

WHITE 


BIAS 

0.15 M RAD 

BIAS 


AHITUDE REFERENCE: 
RANDOM 

0.15 M RAD 

50 MIN 


BIAS 

0.15 M RAD 

BIAS 


ALTIMETER: 

8 METER 

WHITE 

II. 

ENVimriii. 
6E0ID HEIGHT: 

15 METER 

WHITE 


GRAVITY ERROR: 

10 MICRO-G 

2.5 MIN 


DRAG COEFFICIENT: 

10 PERCENT 

BIAS 


EXOSPHERIC TEMP. : 

200 DEGREE K 

50 MIN 




Figure 9. Stellar Update Processing Block Diagram 



Figure 10. Attitude Reference Errors Obtained 
from Covariance Analysis 
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NOTES 

AND S ARE BRIGHT ANO ISOIATEO STARS 

• I OEG^EC RATE (I () lARGE ANCLE NANEUVER 

•360 SEC SEC RATE (I OTRANSIT ACQUISITION 

• TOTAL NANEUVER TINE = 6 * A 2 I ^ • 720 SEC - 12 NIN 

Figure 11, Star Acquisition Attitude Maneuver Profile 


Table 3. Nominal Values of Attitude 
Reference Error Sources 


ERROR SOURCE 

STANDARD DEVIATION 

CORREUTION 

GYRO; 

BIAS 

0.1 DEG/HR 

BIAS 

SCALE FACTOR 

50 PPM 

BIAS 

MISALIGNMENT 

0.05 M RAD 

BIAS 

RANDOM 

0.2 DEG/HR 

WHITE 

STAR SENSOR: 

BIAS 

0.05 M RAD 

BIAS 

RANDOM 

0.05 M RAD 

WHITE 
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THE OPERATIONAL FEASIBILITY OF ORBIT AND ATTITUDE 
DETERMINATION FOR THE GEOSTATIONARY OPERATIONAL 
ENVIRONMENTAL SATELLITE (SMS/GOES) 

USING ONLY IMAGERY DATA 

E. Mack and M. Jurotlch 
Presented by B. Remondi 


U.S. Department of Commerce, NOAA 
National Environmental Satellite Service 


ABSTRACT 

The National Oceanic and Atmospheric Administration (NOAA) uses 
Trilateration Range and Range Rate (TRRR) data, infrared (IR) earth- 
edge data and landmarks for the determination of orbit and attitude 
in the SMS/60ES operations. For many reasons, NOAA would like to 
remove its dependence on TRRR data and determine the orbit and 
attitude using only imagery data. Consequently, NOAA has undertaken 
an investigation to determine the operational feasibility of 
determining orbit and attitude for the SMS/GOES spacecraft using 
imagery alone (either with landmarks only or with the combination of 
landmarks and IR earth-edge data). There are three aspects to this 
investigation: (1) determining the orbit/attitude state under normal 

(no maneuvers) situation, (2) determining the orbit/attitude state 
after the maneuver, and (3) determining the criticality of both quality 
and distribution of the landmark data. 


I . INTRODUCTION 

NOAA currently uses TRRR data, IR earth-edge data, and landmark data 
extracted from visible earth images generated by the on-board Visible 
and Infrared Spin-Scan Radiometer (VISSR) for the determination of 
orbit and attitude in the SMS/GOES operations, NOAA would like to 
remove its dependence on TRRR data and determine the orbit and attitude 
using only imagery data (landmarks only or landmarks with earth-edge 
data) for several reasons: (1) the avoidance of costly processing of 

the TRRR data type and (2) orbit and attitude state could be considered 
a by-product of already existing hardware/software systems. Consequently, 
NOAA has undertaken an investigation to determine the operational 
feasibility of determining orbit and attitude for the SMS/GOES space- 
craft using imagery alone. Three aspects of this investigation are 
as follows: (1) determininq the orbit/attitude state from imagery only 

under normal (non-maneuvers) situation, (2) determining the orbit/ 
attitude state from imagery only after c nunuever, aiio (3) determining 
the criticality of both quantity and distribution of the landmark data. 
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II. LANDMARK AND EARTH-EDGE EXTRACTION AND 
OPERATIONAL DATA PROCESSING PROCEDURES 

There are many different resolutions of SMS/GOES data available. The 
imagery currently used in the orbit and attitude (0/A) determination 
consists of 4 km X 4 km visible and 8 km x 4 km IR. As aprt of the 
ground processing, the earth edges are ascertained using thresholding 
logic and then stored in the IR documentation. (There are 130 bytes 
of documentation data attached to each record of IR data.) They 
are simply the scan (J) and element (I) of the elements of IR data 
at the boundary between earth and space. 

At the end of the ground processing system are ingest computers which 
store these data (imagery plus documentation). The data are moved 
onto the VISSR Data Base (VDB). During this process, the IR earth- 
edge data are extracted. Later, the 0/A models will access these 
data. This data base includes 12 complete visible images of the earth 
and spans approximately six hours centered at spacecraft noon. 

The NOAA Man-Machine Interactive Processing System (MMIPS) access this 
data base to retrieve one picture at a time. These pictures are dis- 
played on a screen and the light pin is placed on a recognizable land- 
mark feature. The scan and sample of this landmark are printed 
automatically. This process is repeated throughout this and other 
visible imagery frames until sufficient landmark data are available 
for 0/A determination. The last step in this process is to add time 
and beta information to these landmarks (I,J) pairs. Thus, an earth- 
edge file and a landmark file have been established from which to 
determine the 0/A state. 

In NOAA's present operation, these imagery data are com; miented with 
TRRR data. This consists of simultaneous ranging data ftjm three 
ground sites. Two of these are unmanned and remote; one of them is 
the prime ground station at Wallops, Virginia. 

Presently, TRRR and landmark data are input into the NOAA 0/A model 
(6E0DYN). Once the orbit and attitude are recovered, the orbit 
(not the attitude) is input along with the ?arth-edge data into the 
the attitude model (PICATT) and an improved attitude solution is 
obtained. The NOAA SMS/GOES ephemerides are then generated using 
these solutions. Nearly all user requirements for ephemeris data are 
satisfied using these ephemerides. 
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III. EXPERIMENTAL RESULTS 

The following results were obtained from three separate spacecraft 
(GOES-1, GOES-2 and SMS-2). The GOES-2 spacecraft replaced the GOES-1 
spacecraft as the operational east spacecraft on August 15, 1977. SMS-2 
is the operational west spacecraft. The coverage from these two space- 
crafts is shown in Figure III.l. 

A. Determination of the Orbit and Attitude State Under Normal (Non- 

Maneuver) Situations 

The GOES-1 and GOES-2 spacecraft were used in this phase of the 
investigation. The subsatellite position for both of these spacecraft, 
at the time of the investigation, was approximately 75° West. The data 
span covered for GOES-1 is from June 23, 1977 thru July 31, 1977. The 
data span covered for WlES-2 is September 18, 1977 thru September 21, 

1977. From Figure III.l, we can see that there exists a number of well- 
distributed lanctaiarks from the imagery data for these spacecraft. 

GQES-1 - The procedure used in this phase of the investigation was 
as follows: starting with an operational a priori estimate and seven- 

days .vorth of landmarks from the Man-Machine Interactive Processing 
System (MMIPS), the GEODYN orbit determination system was used to deter- 
mine both orbit and attitude from landmarks-only. Next, using the 
landmark-only orbit and the IR earth edge data in the Horizon Picture 
Attitude Program (PICATT), we compute a second attitude solution. Thus 
in reality, we have two daily attitude determinations where the orbits 
for both solutions are the same. The two solutions described above are 
then used in the Gridding Error Assessment System to compute the x- 
direction shift (east-west shift), y-direction shift (north-south shift), 
and rotation which are used to judge the quality of the solution. The 
first two of these give an indication of what the grid error would be 
is these solutions were used in the gridding. 

All the landmark-only solutions in this phase of the investigation were 
determined from seven-days worth of landmarks. The resultant orbit, 
attitude and grid errors were compared with the operational orbit, 
attitude, and grid errors and the following results were obtained: (a) 

Figure III. 2 and Figure III. 3 show the east-west and north-south grid 
error produced from each solution in terms of grid errors in the landmark- 
only orbit solution with the PICATT attitude solutiw. Table III.l shows 
the mean and standard deviations for the grid errors computer by the 
Gridding Error Assessment System for the three solutions described above, 
(b) In addition to looking at the grid error produced from each of the 
solutions, we can also examine the orbit and attitude solutions to see 
what differences existed. Table III. 2 shows the mean orbital differences 
based upon the sixteen orbit solutions used to produce the grid error in 


- 
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Figure III. 2 and Figure III. 3. Tables III. 3 shows the mean attituae 
differences based upon the sixteen attitude solutions used to produce 
the grid errors in Figure til. 2 and Figure III. 3. From Table III. 3, we 
can see that the attitude solution is considerably improved by the 
addition of the IR earth-edge data in the attitude solution. 

TABLE III.l 


GRiID ERROR SOURCE 

X-OIRECTION SHIFT 

Y-DIRECTION SHIFT 


Mean 

Std. Dev. 

Mean 

Std. Dev. 

Landmark Only 
Orbit & Attitude 

-5.724 km 

3.537 km 

-15.2357 km 

5.9745 km 

Landmark Only 
Orbit & PICATT 
Attitude 

-5.2097 km 

4.1820 km 

-12.1983 km 

3.3229 km 

Operation 
Orbit & Attitude 

-5.3357 Ian 

7.6966 km 

-10.6403 km 

6.2384 km 


TABLE III. 2 
Orbital Differences 

Landmark Only vs. Operational Solutions 


Semi -Major Axis 
Eccentricity 
Inclination 
Mean Longitude 


Mean 

233.25 meters 
0.00001 
0.008 deg 
0.008 deg 


Standard Deviation 
5.00 meters 
0.00003 
0.005 deg 
0.007 deg 


TABLE I I I. 3 
Attitude Differences 

Landmark Only vs. 
Operational 

Mean Std. Dev. 

3.86 deg. 6.65 deg. 

.08 deg. .01 deg. 


Landmark and Earth- 
Edge vs. Operational 

Mean Std. Dev. 

2.62 deg. 

.02 deg. 


Spin Axis Right Ascension 
Spin Axis Declination 
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7.78 deg. 
.03 deg. 



2 . GOES-2 - On September 16, 1977, a landmark-only orbit solution with 
PICATT attitude was computed and used in the computation of the opera- 
tional grids used by Wallops CDA to grid the images fran GOES-2. The 
operational grids are computed in three eight-hour periods as follows. 

gridding period 1 = 0230Z thru 1030Z 

gridding period 2 = 10302 thru 18302 

gridding period 3 = 18302 thru 02302 

Moreover, only one period of grids is computed per day so that the 
lanctaark-only orbit computed on September 16, 1977, was used in the 
following operational grids: 

Septanber 18 for 02302 thru 10302 
September 19 for 02302 thru 18302 
September 20 for full 24 hours 

A number of pictures gridded during those periods were examined and the 
grid error was measured by hand. Figure I I I. 4 and Figure I II. 5 shows 
the hand-measured grid error from those periods. It should be noted that 
the grid error can be measured by hand to an accuracy of about 7.5 
kilometers for GOES-2. 

B. Determinaticn of the Orbit and Attitude State After a Maneuver 

The Synchronous Meteorological Satellite (SMS-2) was used in this 
phase of the investigation. The subsatellite position for this space- 
craft is approximately 135° West. The data span covered for SMS-2 is 
June 20, 1977 to June 27, 1977. From Figure III.l, we can see that we 
do not have as many well -distributed landmarks available from the imagery 
data for this spacecraft as we had with the GOES spacecraft. 

The procedure used in this part of the investigation was as follows: 
Starting with the post-maneuver predictions for orbit and attitude on 
June 20 and one day's accumulation of landmarks from the MMIPS, we used 
the GEODYN orbit determination system to determine both orbit and atti- 
tude from landmarks-only. If GEODYN is unable to get a solution, add 
another day's worth of landmarks and repeat the GEODYN run. Continue 
this process each time adding another day's worth of landmarks until we 
are able to get a GEODYN solution. 

Using the above procedure, a successful GEODYN solution was achieved 
with seven day's worth of landmarks. The GEODYN solutions were compared 
with the corresponding operational solutions and the results are shown 
in Table III. 4. 



We can see from Table I I I. 4 that a landmark-only solution was obtained 
with seven days worth of landmark that agrees very closely with the 
operational solutions. This was significant, since with the SMS-2 space 
craft, we do not have well-distributed landmarks (most landmarks are 
along the west coast of the United States). 

TABLE I II. 4 


Landmark Only 
Orbit & Attitude 


Operational Orbit & 
Attitude Solution 


Semi -Major Axis 
Eccentricity 
Inclination 
Mean Longitude 
Right Ascension 
Declination 
Data Span 


42166224.8378m 
0.0001589 
0.044409 deg. 
383.4808567 deg. 
19.869378 deg. 
-89.821767 deg. 

7 days 


42166009.0m 
0.0001524 
0.045367 deg. 
383.560237 deg. 
19.80943 deg. 
-89.84232 deg. 

2 days 


Table III. 5 shows the orbit and attitude differences between the post- 
maneuver landmark-only and post-maneuver operational orbit. 

TABLE III. 5 


Semi -Major Axis 

Eccentricity 

Inclination 

Mean Longitude 

Spin Axis Right Ascension 

Spin Axis Declination 


Landmark-Only vs. Operational 

216.83 meters 
0,000006 
0.009 deg. 

0.08 deg, 

0.06 deg. 

0,02 deg. 


TABLE I II. 6 


Semi -Major Axis 
Eccentricity 
Inclination 
Mean Longitude 
Right Ascension 


Landmark Only 
Solution With 
3 Sets of Landmark 2 

42166224.8378m 
0.001589 i 

0.044409 deg, ' 

383.4808567 deg. 
19.869378 deg. 


Landmark Only 
Solution With 
! Sets of Landmark 

42166389,324m 
0.0001596 
0.058101 deg. 
383.468963 deg. 
19.941449 deg. 


Landmark Only 
Solution With 
1 Set of Landmark 

42186470,3063m 
0.00064219 
0.04461 deg. 
383.515453 deg. 
19.953716 deg. 
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TABLE I I I. 6 CONT. 


Oeci inatlon 
Data Span 

# Elements 

# Lines 


Landmark Only 
Solution With 
3 Sets of Landmark 


Landmark Only 
Solution With 
2 Sets of Landmark 


-89.821767 deg. 
7 days 
126 
126 


-89.849354 deg. 
7 days 
44 
44 


Landmark Only 
Solution With 
1 Set of Lanchnark 

-89.824916 deg. 
7 days 
21 
21 


C. Critcality of Both Quantity and Distribution of the Landmark Data 


With both the GOES and SMS-2 spacecrafts, we try to extract landmarks 
so that we have three separate sets of landmark locations in order to 
enhance the geometry of the solutions. However, only the SMS-2 space- 
craft was used in this phase of the investigation. The procedure used 
in this investigation was as follows; Start with the post-maneuver 
predictions for orbit and attitude on June 20 and seven days worth of 
landmarks with one set of landmarks' locations removed from the land- 
mark data. Then use the 6E0DYN orbit determination system to determine 
both orbit and attitude from landmarks-only. Next, remove a second set 
of landmark locations from the landmark data and compare a second 
GEOOYN orbit and attitude solution. Table III. 6 shows the results of 
these 6E0DYN solutions as compared with the landmark-only solution 
from Section B. 


IV. CONCLUSIONS 

Experimental results from the three data evaluation periods (June 20, 
June 23, Septen6er 18) on the three geostationary spacecrafts (SMS-2, 
GOES-1, GOES-2) have demonstrated that; 

1. Using existing landmark extraction and identification techniques for 
the east geostationary spacecraft, we can maintain a high quality orbit 
and attitude state with imagery data only. 

2. Using existing landmark extraction and identification techniques for 
the west geostationary spacecraft, we can recover a high quality orbit 
and attitude state with imagery data only in approximately seven days. 
Even though this recovery period appears prohibitively large, we should 
note that SMS-2 does not have well-distributed landmarks and the addi- 
tion of IR earth-edge data in the orbit solution should improve the 
recovery time. 

3. Analysis of the landmark data for SMS-2 indicates that the removal 
of two of the three sets of landmarks from the data would seriously 
endanger the ability to achieve stable orbital elements. However, it 
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should be noted that for SMS-2 most of the landmartts extracted are along 
the west coast of the United States and once the set of west coast land- 
marks are removed, the geometry is seriously hampered. 



I'igurt* III.l. GOES Operational Coverage. 









attitude 



Figure m. 3. Norlh-South Grid Error Comparison of Orbit and Attitude 
Deter minati<m Methods. 
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RECURSIVE ESTIMATOR FOR OSO-8 
ATTITUDE 

Robert D, Headrick and Duke Y. Park* 

Computer Sciences Corporation 

INTRODUCTION 

Hie validation and early results of the Recursive Estimation Attitude Program 
(REAP) were presented in an earlier paper (References 1 and 2) in the May 1976 
Estimation Theory Symposium. This paper presents modifications and en- 
hancements that have been made to the program since then. 

The REAP program is used for research and special production for definitive 
attitudes for Orbiting Sol^r Observatory (OSO)-8. The objective is to deter- 
mine continuous attitudes to ± 0. 05 degree accuracy from Sun and star slit sen- 
sors mounted on the spinning portion of the spacecraft. The bulk of the attitude 
production is performed by a Weighted Least Squares (WLS) batch processor, 
but REAP is used for problem passes such as those involving gas jet maneuvers, 
sparse star fields, or star sensor saturation by high energy particles in the 
South Atlantic Anomaly. 

The star sensor has a near-vertical azimuth slit and a canted elevation slit, 
with a single photomultiplier tube as a detector. The Sun sensor similarly has 
a vertical and a canted slit. Thus, the sensings are in the form of time-tagged 
events, where the time is taken from the spacecraft clock pulses. Other sen- 
sors of lower accuracy are also available: a single-axis magnetometer which 
given the time of the rising and falling nulls and a Shaft Angle Encoder (SAE) 
which gives the relative azimuth between the sail (held approximately normal 
to the sunline) and the wheel (rotating at about 6 rpm). 

♦Work supported by Goddard Space Flight Center, National Aeronautics and Space 
Administration under Contracts NAS 5-11999 and NAS 5-24300. 
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A diagram of the spacecraft is given in F^ure 1. The attitude angles are re- 
ferenced to the Geocentric Solar Ecliptic (GSE) coordinate system in which the 
axis is pointed toward the center of the Sun, z_ is to the ecliptic north 
pole, and y completes the right-handed system. The wheel (or body) 

axis orientation is described in terms of roll and pitch angles from these 
reference axes* The roll angle, 0 , shown negative for clarity. Is taken about 



Xg,Yg,Zg - ECUPTIC COORDINATES 
WHEEL COORDINATES 


A - ROLL ANGLE 

n - PITCH ANGLE 

$ ^ ASPECT ANGLE 


Figure 1. OSO-8 Orientation Angles 



the sunline* and the pitch angle, n » is about the intermediate axis, y^. The 
final rotaticHi is about the body axis through the aspect angle 0. 

Mathematical Formulation 

The computational sequence, which is given in Figure 2, is the same as in the 
earlier versions. The mathematical formulation, however, has been simplified 
considerably from the original state representation which included both the 
attitude quaternion and the angular momentum vector. Since there is no evi- 
dence of nutation, this representation was clearly redundant. Now the state 
vector, X . is given by 

X= »» 

-W. 

where ^ = roll angle 
7J= pitch angle 
fi- aspect angle 
w- spin rate 

The measurement error is given by 

Ay^ = Oj • = 0^ D W 

- th - . . 

where U is the slit normal vector for the j event, V is the reference vec- 
tor transformed to the body system, W is the reference vector ’n the ecliptic 
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Figure 2. Sequential Estimator Computational Sequence 


















(GSE) system, and D is the attitude direction cosine matrix referenced to the 
GSE system, given by 


/ C^CT) 

cj3st)s^ + s/3c0 

-c)5st)c0 + 

S^S0 

j-SjScT) 

-S/SsT)S0 + CjSc0 

S^ST)C0 + 

C$S0 

\ St) 

-Ct)s0 

Ct)C0 



In the predictive step the attitude dynamics have been simplified to the spin-axis 
approximation, where it is assumed that the body z-axis is aligned with the 
angular momentum vector (i.e. , no nutation occurs). The angular momentum 
vector L. is computed at each step for use in predictions, but it is not included 
In the state estimation procedure. It is updated by the equation 

L. , =L. + Aty^Torques. - m XL. 
j+1 j * j ‘Sun j 


where is in GSE coordinates. The various torque models, which are av- 
eraged over a spin period, are now computed in the GSE system to avoid con- 
versions. The torque models include 


• Gravity-gradient 

• Solar radiation 

• Magnetic residual and control torques 

• Attitude control jets 


The uL term accounts for the rotation of the GSE coordinate system. The 
Sun 

state vector is predicted by 


j+1 




where v, 


J+1 


is a random noise vector. 
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The state transition matrix, 0, is given by 




&x. 


J 


0 

1 

0 

0 


0 

0 

1 

0 


0 

0 

At 

0 


The dtynamic noise term is given by 


Q = 


0 

0 

0 


0 

0 

0 


0 

0 

Q, 


0 

0 

0 

Q, 


where and are constants which account for the immodeled torques and 
computational noise. 

The measurement sensitivity matrix, H, is given by 


« _ _tj^3DW 

n - ^ A * 


bx 


J 


asc 


where the attitude dependence is contained solely in the direction cosine matrbc. 
To illustrate the simplicity of these partlals, we develop them explicitly as 


d0 


0 C$snc0~8ps<f> C/5STJS0+ s^ctf 

= 1 0 -s)3sT)c0-ci3s0 -s^sT)s0 + c;5c0 


-CT|C0 


-CT)80 
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3D 

an 


-0/5 sn 

c0crfa0 

c$crfc0 

s08ri 

-s/Scns0 

80cric0 

cn 

81)80 

-snc0 


/ 


IS 

a/5 



-Sj3sT}S0+ Cj8c0 
•~C$8l)S0” S$c0 
0 


s/3snc0+ c^s0 
C0srjc0 - 8030 
0 


The U and W vectors which pre-and post-multiply the partials SD/Sx serve 

to select and weight the terms in the mafrbc. For example, if we have a Sun-slit 

— T 

in the body x-z plane, the slit normal is U = (0, 1, 0) . The Sun vector in 
GSE coordinates is 



Thus, when we carry out the multiplication, we select the (2, 1) term of the 
SD/ax matrices to give 


H = ^0, S0SP. - C0CT), Oj 

These measurement sensitivity matrices are simple enough to allow checking 
by hand computations. 
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Note that the roll partial is identically zero since the vector product 
(dD/d0) W is zero. Thus, there can be no direct update of the roll angle on 
Sun dtta, no matter where the Sun slit is located. Since the a priori covariance 
matrix is diagonal, the initial roll gain is also zero. When the data span 
begins in Sunlight, the roll update remains zero throughout orbit day. After 
the star data comes in at orbit night, however, the off-diagonal correlation 
terms in the covariance matrix build up. During the next orbit day these give 

rise to a nonzero K . and provide an indirect update to the measured roll angle. 

w 

Data Editing 

Ibe method of star Identification is the direct match technique, which selects 
out of the catalog that star which has the smallest residual errors as the ref- 
erence star. In order to avoid misidentlfications from the many false trigger- 
ings which occur, considerable care is taken in data editing. The total star 
catalc^ is reduced to stars brighter than 3. 5 magnitude In the swath swept out 
by the star sensor field of view about the a priori spin axis. Events are edited 
out for the following conditions: 

• Events outside nominal azimuth (counts do not agree with the 
spacecraft clock) 

• Events outside the elevation limits 

• Triple-crossing flag in telemetry indicates 3 events (and possible 
ambiguity) within 125 milliseconds. 

• Sensor saturation due to the South Atlantic Anomaly, indicated by 
more than 128 events during a major frame ot 20. 48 seconds 

• Stars occulted by the Earth 

• Stars which are not identified as occurring from both the azimuth 
and elevation slits. Since the events are considered singly, special 
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contingent logic Is required to hold all candidate slit 1 reference 
stars and states until the data from slit 2 is processed. 

* Events outside fixed tolerance limits 


• Events outside adaptive limits. This test has been added since the 
the last paper to provide T) -sigma data rejection based on the stan- 
dard deviation of the residual 


0 =v^PH^ + R 
AY ' 


If Ay >na^y t the data is rejected, where n is an input parameter. 


Results 


Results of the new 4-component model are shown in Figure 3 compared to the 
earlier solutions, labeled V.2 and V.4. The spin rate agreed so closely with 
the earlier results that it was not plotted. The oscillations in spin rate at the 
orbit period are due to thermal expansion. As the spacecraft enters sunlight, 
it expands, increasing the moment of inertia. Since angular momentum is 
conserved, the spin rate must cecrease. The opposite effect occurs when the 
spacecraft enters orbit night. 

The roll and pitch angle results are shown, where the horizontal bars indicate 
the allowed error. It can be seen that all versions agree to within ± 0. 03 de- 
gree. The 4 -component pitch solutions show a small systematic oscillation at 
the orbit period, which was not predicted by the torque models. When the data 
is turned off, the predicted state follows the dashed (V.2) curve very closely. 
The cycles at twice the orbit frequency are characteristic of gravity-gradient 
torques. The new observed pattern is indicative of aerodynamic torques from 
the fixed sail. These torques were not modeled originally because they were 
known to average out over an orbit period (Reference 3) and thus not lead to any 
secular error growth. Also, the aerodynamic terms are very complicated to 
model accurately. 
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Figure 3. Attitude Results on August 12, 1975, Compared with 
Earlier Versions and Weighted Least Squares Results 
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To test whether these oscillations were due to aerodynamics, we incorporated 
a crude aero torque model into the prediction. This model includes only the 
sail and assumes a constant density over the circular orbit. The predicted 
state (with aero torques but without data) agreed very well with the estimated 
state solutions (with data but without aero torques). It shows the same pattern 
of alternating deep and narrow valleys, though the deep valleys were not as 
deep as in the estimated solution. 

Conclusions 

The oscillations at the orbit period have been shown to be real, due to un- 
modeled aerodynamic torques. These torques may slightly degrade the accuracy 
of the solutions during passes with minimal data, but are not required to be 
modeled for data spans with dense data as in Figure 3. Most of the data proc- 
essed so far has been with Version 7, which agreed with die 4-component solu- 
tions to within 0.02 degrees though it was relatively closer to the V.2 curve. 

This solution was not plotted to avoid over-complicating the graph. 

Equally important, the sensitivity of the 4-component version to unmodeled 
torques shows that the filter is 'open* to new data. The main parameters to be 
tuned were the state noise terms and , which were adjusted to make 
the attitude standard deviation from the covariu je matrix approximately equal 
to the average deviation of the residuals. It further shows that the earlier ver- 
sion (V. 2) was closed and did not follow the data. Version 4 was more open ^ 
but the solution also showed artifacts, such as those at 1700, which were duu 
to the multiple path logic. 

bi addition to being easier to tune than the earlier versions, the 4-component 
model is also twice as fast. It runs at a ratio of 50:1 of real-time to CPU 
time on the Univac 1108. 
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PERFORMANCE OF GROUND ATTITUDE 
DETERMINATION PROCEDURES 
FOR HE AO-1* 

Lawrence Fallon III and Conrad R. Sturdi 

Computer Sciences Corporation 
Silver Spring, Maryland 


INTRODUCTION 

HEAO-1 was placed in a low circular orbit on August 12, 1977. It wei^s 
approximately 3150 kilograms (7000 pounds) and is the heaviest unmanned 
Earth-orbiting satellite launched by NASA to date. The observatory will sur- 
vey and map X-ray sources throughout the celestial sphere and will measure 
low energy gamma-ray flux. HEAO-1 is controlled to scan about the sunline 
and thus will provide a complete survey of the sky in 6 months. Ground attitude 
support is provided at GSFC by the HEAO-1 Attitude Ground Support System 
(AGSS). Information telemetered from Sun sensors, gyroscopes, star trackers, 
and an onboard computer are used by the AGSS to compute updates to the on- 
board attitude reference and gyro calibration parameters. The onboard com- 
puter utilizes these updates in providing continuous attitudes (accurate to 
0. 25-degree) for use in the observatory's attitude control procedures. This 
paper will discuss the relationship between HEAO-1 onboard and ground proc- 
essing, the procedures used by the AGSS in computing attitude and gyro cali- 
bration updates, and the performance of these procedures in the HEAO-1 
postlaunch environment. 


♦Work supported by the Spacecraft Control Programs Branch, Goddard Space 
Fli^t Center, National Aeronautics and Space Administration, imder Contract 
NAS 5-11999. 
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Figure 1. HEAO-1 Onboard and Ground Processing 


HEAO-1 ONBOARD AND GROUND PROCESSING 

The prime contractor of the HE AO mission, TRW Systems Group of Redondo 
Beach, California, devised the onboard system of attitude determination and 
control shown schematically in Figure 1 (Reference 1). This system consists 
of four Bendix rate-integrating gyroscopes, coarse and fine Sun sensors manu- 
factured by TRW, two star trackers manufactured by the Ball Brothers Research 
Corporation, a Control Data Corporation (CDC) flight computer, and reaction 
control thrusters. The gyros measure angular displacements of the observa- 
ory during each 320-millisecond minor frame cycle. The input axes of the 
four gyros are skewed with respect to each other so that complete 3-axis ro- 
tational information is provided by any three of the four gyros. Output signals 
from the gyros are received by the gyro processing module in the fli^t com- 
puter. In this module, a nominal scale factor (K) and a 3-by-4 scale factor 




correction and alignment matrix [r 1 ui-e implied to the gyro outputs to com- 
pute an angular velocity in the spacecraft reference fitune. This angular ve- 
locity is further corrected for gyro null shift using a drift rate vector (S). 

The resultant angular velocity (to ) is available each 320-millisecond minor 
frame to the control law and to the attitude propagation module. 

An incremental rotation quaternion (corresponding to the observatory's rotation 
during the 320-millisecond minor frame) is constructed in the attitude propa- 
gation module using t* e angular velocity. The observatory's estimate of its 
attitu(te is then propagated through the minor frame via quaternion multiplication. 
This propagated attitude (Q) is then available every minor frame to the control 
law, and every 40. 96 seconds in telemetry. 

In the celestial scan control law, the angulpr velocities and propagated attitudes 
are compared to a commanded scan rate (u)^) and a target attitude (Q^) to 
generate thruster activation commands for scan rate and Z-axis attitude con- 
trol. (A celestial point mode to be entered occasionally later in the mission 
compares Q and to generate activation commands for 3-axis attitude 
stabilization. ) In the celestial scan mode, target attitudes are sent to the ob- 
servatory twice a day to cause its Z-axis to follow the sunline within 1 degree. 

Errors in the gyro calibration parameters, [t 1 and b, cause errors in the 
angular velocities and hence errors in the propagated attitudes. To prevent the 
divergence of the observatory's Z-axis from the Sun and to initialize the attitude 
reference following launch, an attitude reference update procedure was ^vised 
by TRW. In this procedure, corrections are applied to the attitude reference 
based upon ground attitude solutions computed by the AC^ using the telemetered 
star tracker data. Corrections may also be applied using the two-axis angular 
displacements of the Sun from the observatory's Z-axis, as measured by the 
fine Sun sensor. This type of correction proviues no attitude improvement about 
the sunline, but is sufficient to cause the observatory' ^ Z-axis to follow the 
Sun. To further improve the accuracy of the onboard attitude reference, the 
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AGSS periodically transmits new values for the gyro drift rate vector, b. The 
AGSS estimates new drift rates by monitoring the divergence of the propagated 
attitudes from the attitude solutions computed using the star tracker data. The 
capability for correcting the calibration matrix, [r1 , also exists; but as long 
as 5 is nearly constant, corrections to both b and [t] are not simultaneously 
observable. For this reason, it is unlikely that new values for [t ] will be 
sent to the observatory until later in tlie mission. In the attitude control en- 
vironment, the AGSS is required to compute reference and gyro calibration 
updates of sufficient accuracy and frequency that the observatory's attitude 
reference is maintained within a 0. 25-degree accuracy level. 

FUNCTIONAL OVERVIEW OF THE AGSS 

The AGSS was devised by GSFC an< Computer Sciences Corporation (CSC) to 
meet the HEAO-1 attitude support requirements. A functional description of 
the AGSS is presented in Figure 2. Telemetry from gyros. Sun sensors, star 
trackers, and the onboard computer is processed, and necessary star catalogs 
are generated. Stars are selected from a SKYMAP run catalog to form a sub- 
catalog which is typically a 12-degree-wide band orthogonal to the Sun vector 
(Reference 2). The major attitude processing sequence of the AGSS is then in- 
voked. The first step in this procedure is to define inteiwals based upon data 
continuity and quality. Minor frame data is tlien processed. In this step, atti- 
tudes are propagated to desired times using the gyro data. Star tracker data 
is edited and calibrated to form observed star unit vectors in the spacecraft 
reference frame. These vectors are then transformed into a common reference 
frame using the propagated attitudes at the time of each star tracker observa- 
tion. Average unit vectors are then computed using approximately 20 sequen- 
tial sightings of each different star. WTien this procedure has been completed 
for all star tracker observ'atlons in a particular interval, attitude solution proc- 
essing is initiated. If an adequate attitude estimate is not available, an attitude 
acquisition procedure, described below, is used to provide a 3-axis attitude. 






When a suitable estimate has been obtained, small segments or "snapshots** of 
star observations are identified and least-squares attitude solutions are calcu- 
lated using a procedure also described below. Following computation of several 
attitude solutions spaced in time throughout the interval of data, an appropriate 
solution is selected to form an attitude reference update. Solutions and associ- 
ated propagated attitudes are then assembled and stored for later use in gyro 
calibration. After this attitude processing sequence has been conq>leted for 
approximately one orbit of data, a short term drift rate solution is calculated. 
After several orbits have been processed, the ^rro calibration module (see be- 
low) is again invoked to estimate a longer term drift rate solution for possible 
transmission to the observatory. The calibration matrix, [r1 , could also be 
estimated at this time if sufficient observability is anticipated. 

A. Initial Attitude Acquisition 

The major inputs for initial attitude acquisition are observed Sun and star unit 
vectors transformed to a common reference frame and a star catalog in the 
shape of a band orthogonal to a Sun vector, supplied by ephemeris. This catalog 
has a limiting magnitude 0. 7-magnitude fainter than the star tracker threshold; 
the stars are arranged in order of azimuth about the catalog pole. Star identi- 
fication begins with the selection of a key observation, nominally the ftrst 
observation with intensity brighter than an input limit. The key obseirvation 
is matched to each observation in turn, yielding tentative X-axis phases. In 
each case candidates for the other observations are formed by matching the 
azimuthal difference between observation and ca:alog stars with the azimuthal 
difference between the key observation and the star tentatively associated with 
it. Angular separations between the key observation and the other observations 
are compared to the separations between the star matched with the key obser- 
vation and the above candidates. The candidate ^f any) which best satisfies this 
angular separation within a small tolerance is identified with the observation. 
These identifications are verified by demanding that the angular separation 


224 



between the observation and the observed Sun vector match the separation 
between the catalog star and ephemeris Sun vector. The largest such set of 
identifications is considered to be correct. A least-squares attitude solution 
is then compited with this set of identifications using a batch least squares 
algorithm developed by P. Davenport (Reference 3). To assess the quality of 
this 3 -axis attitude solution, the procedure is repeated with different key ob- 
servations. 

B. Snapshot Attitude Determin ation 

The initial attitude estimate is used to transform the observed star unit vectors 
to an approximate Geocentric Inertial (GCl) reference frame. A short segment 
or "snapshot" of the transformed observations is then selected for determina- 
tion of an improved attitude estimate in the form of a correction quaternion. 

The catalog is searched over small intervals of azimuth centered on each ob- 
servation for stars whose angular separation from the observation is less than 
an input tolerance. A noncolinear triplet of observations, with separations 
exceeding specified limits and a minimal number of candidate stars, is selected. 
Identification of this triplet is made by matching the angular separations of the 
three observations with corresponding separations of the candidate stars. The 
other snapshot observations are identified with the candidates whoso angular 
separations from the identified triplet stars best match tlie corresponding 
separations between these observations and the triplet of observations. Finally, 
a least-squares attitude solution for the mid-snapshot time is obtained from the 
weighted values of the observation unit vectors and the unit vectors of the 
corresponding catalog stars. 

C. Gyro Calibration 

Input to the gyro calibration module consists of pairs of sn^shot attitude solu- 
tions (Qgx> Qs2^j associated propagated attitudes (Qpi* Qp2)j the two 
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attitudes in each pair are typically separated in time by approximately 30 to 
90 minutes. Also provided are continuous histories of propagated attitudes 
selected every 41 seconds between each (Qp]^)^ and (Qp 2 )j values of 

the gyro calibration parameters, [r] and b , which were used in calculating 
each set of propagated attitudes. 


Several sn^shot pairs are selected, provided that the same [r] and b values 
were used in creating the associated propagated attitudes. At least one pair 
of attitudes is necessary to estimate b and at least four pairs are needed to 
estimate both b and [t]. Corrections to b and, optionally, [r] are com- 
puted using a batch least-squares algorithm developed by P. Davenport 
(Reference 4). In this procedure, the following loss function is minimized: 


J= £ Cz -[H]X]'^[Z -Ch]x] 
j = 1 J •» J J 


where N is the number of snapshot pairs selected; the Z. are the vector 

-1 ^ -1 - 

parts of the error quaternions given by Qg^ Qgg Qp^ ); X is a 
12-vector containing corrections to b and [t]; and the [H ] are 3 by 12 
matrices containing the partial derivatives of the with respect to X. The 
[H.] are calculated sequentially using the propagated attitude histories. The 
Z^ would be zero if no errors existed in [t ] or b. 

AGSS PERFORMANCE FOLLOWING LAUNCH 


A. Initial Attitude Acquisition 

The hi^ voltage was applied to the star trackers two days following launch. 
Initial attitude acquisition was obtained from the data described in the second 
column of Table 1. The prelaunch gyro calibration parameters were used 
during processing of the data. The accuracy of the 3-axis attitude solution, 
estimated to be dsO. 2-degree per axis, was sufficient for subsequent snapshot 
processing. 
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Later in the day» erroneous gyro calibration parameters were transmitted to 
the spacecraft. The use of these parameters by the onboard computer resulted 
in a large error in the propagated attitude and necessitated a second initial 
attitude lusquisition attempt described in the third column of Table 1. The 
star tracker direshold had been commanded to the sixth-magnitude level and 
the calibration lights turned on. In this mode, observations of calibration lights 
account for about one-half of the total; these observations are rejected before 
star identification is attempted. Because the X-axis phase was known to within 
±20 degrees, only one-ninth of the band star catalog 7-magnibade-limlt) was 
used for key-observation matches. Computation time was also shortened by 
requiring idenfification of only the first 30 star observations. The higher accu- 
racy of the second attitude solution was due to the use of gyro calibration param- 
eters estimated in flight from data following the first attitude acquisition. 

B. Snapshot Attitude Determination 

Snapshot processing statistics from a typical orbit of data are displayed in 
Table 2. About two-thirds of the original observations were rejected because 
they were caused by calibration li^ts, were excessively noisy, or lay outside 
the usable star tracker field-of-view. Nearly all of the observations for which 
identifications were attempted were actually identified. Most of the identifica- 
tions which were subsequently rejected in the least-squares procedure were 
rejected due to the presence of a neighboring star within a specified angular 
separation range. The RMS of the angular residuals between the transformed 
observations and associated catalog stars is given in the last row of Table 2. 

C. Reference and Gyro Calibration Update Procedures 

The performance of the reference and gyro calibration update procedure is best 
assessed examining the divergence between the onboard attitude reference 
and AGSS sniqishot solutions as a function of time. Figure 3 illustrates the total 
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IMATEO ERROR, EACH AXIS (DEG) ±0.2 ±0.03 
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Figure 3. Attitude Propagation Accuracy - Week Following 8/14/77 






Src difference between AGSS snapshot solutions and onboard attitudes at the 
same time» for the week following star tracker turn-on. .After initial acquisi- 
tion and until the first reference update was received, the onboard attitude was 
in error by approximately 60 degrees (nearly all in phase). When this first 
reference update was sent, an erroneous drift update was also transmitted. 

The attitude reference was initially corrected, but it then began to diverge 
rapidly to an error of nearly 16 degrees (when a correct drift update was re- 
ceived the observatory). The next reference update caused the onboard 
attitude error to decrease within the 0. 25-degree level, Additio*'-’! reference 
and gyro updates further improved onboard accuracy throughout the remainder 
of the week following August 14. 

Onboard propagation accuracy for the second week following star tracker turn- 
on is illustrated in Figure 4. Excluding the large errors observed ^ .August 22 
(caused by an unusually long interval when no data was received by the AGSS), 
additional improvement in onboard propagation accuracy is obseiwed. Onboard 
performance after several weeks closely parallels the first half of the week 
following September 16, as shown in Figure S. In this period, onboard accuracy 
is held to near or better them the 0. OS-degree level (the level r^^ired for 
post-facto definitive attitude determination). On September 21, however, a 
commanded scan rate change caused a rapid increase in onboard error due to 
the strong dependence of drift rate solutions on the scan rate. A new drift rate 
was estimated using data following the scan rate change, and sent to the ohser- 
vatorj' on September 22. Propagation accuracy then returned to within the 
0. 05-degree level. 

The gyro drift rate solutions for the week following September 16 are shown in 
Figure 6. Short-term drift solutions (obtained using data from approximately 
one orbit) are connected by a bold line. The apparent variation of the X and Y 
components is due more likely to difficulties in observability than to true drift 
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Figure 5. Attitude Propagation Accuracy - Week Following 9/16/77 






variation. The values sent to the observatory are illustrated by a thin horizon- 
tal line. These values were estimated using data acquired from several orbits 
in the 24 hours prior to transmission. The time variation of these longer-term 
solutions Is far less dramatic. 

CONCLUSIONS 

The postlaunch performance of the ACiSS may bo summarized as follows: 

• Initial attitude acquisition was achieved with sufficient accuracy for 
subsequent snapshot processing. 

• Snapshot processing provides reliable attitudes for reference up- 
dates, gyro calibrations, and quality assurance of onboard attitude 
propagation. 

• Attitude updates and gyro calibrations provided by the AGSS enable 
the onboard computer to maintain its attitude reference to well 
within the 0. 25-degree accuracy requirement (and often within the 
0.05-degree definitive accuracy requirement). 

The success of the AGSS in its ^xistlaunch supjxjrt is due in part to tlie perform- 
ance of the observatory'* s attitude sensors, particularly tlie gy'ros. The com- 
bined performance of the various components of HEAO-1 onboard and ground 
processing demonstrates the technical feasibility of using an onboard computer 
to supply attitudes of definitive quality. 
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ISEE-C ATTITUDE DETERMINATION USING 
FINE SUN SENSOR DATA ONLY* 


Lawrence P. Gunshol 
Computer Sciences Corporation 


ABSTRACT 

The liitemational Sun-Earth Explorer-C (ISEE-C) will be spin stabilized with 
the spin axis attitude nominally maintained within 1.0 degree of the north eclip- 
tic pole. ISEE-C will be stationed in the vicinity of the Sun-Earth interior 
libration point where the only significant attitude perturbation is due to solar 
radiation pressure. The primary attitude sensors are two fully redundant 
Fine Sun Sensors (FSSs). The specified sensor accuracy after calibration is 
±0.03 degrees for a ±3 degree range about the spin plane. An operational 
requirement is to determine the attitude to within 1.0 degree using up to 30 
days of FSS data. To do this, we have developed techniques to determine the 
spin axis attitude usii^ FSS data only. At any given time, the Sun angle speci- 
fies the orientation of the spin axis relative to the sunline. The instantaneous 
time rate of change of the Sun angle is directly proportional to the orientation 
of the spin axis relative to a reference plane that is normal to the ecliptic. 
Thus, the spin axis attitude can be determined when sufficient data has been 
collected to accurately measure the rate of change of the Sun angle. The un- 
certainties can be computed directly from the uncertainties in the coefficients 
of the smoothed Sun angle curve. 

The FSS-only technique is unique in that ephemeris vectors are required only 
to transform the attitude results to more conventional coordinate frames. The 
combination of the mission geometry and the FSS accuracy make ISEE-C an 
ideal mission for appl 3 dng this method. However, the technique can be used 
on other missions, such as spin stabilized geosynchronous missions. 


* 

This work was performed under Contract NAS 5-11999 with NASA Goddard 
Space Flight Center. 



ISEE Program Overview; The International Sun- Earth Explorer (ISEE) Program 
is an international cooperative program which will use three coordinated space- 
craft to advance knowledge of the magnetosphere, interplanetary space, and 
their interactions. On October 23, 1977, the IBEE-A (mother) and ISEE-B 
(daughter) were launched aboard the same Thor Delta launch vehicle into a 
highly eccentric Earth orbit (apogee equal to 22 Earth radii, perigee equal to 
approximately 300 km). The third spacecraft, ISEE-C, will be launched anoai'd 
a Delta 2914 rocket and will be placed in a heliocentric orbit near the unstable 
libration point between the Earth and the Sun, a distance of approximately 
4 lunar orbit radii or 0.01 astronomical unit from the Earth. The National 
Aeronautics and Space Administration has responsibility for development of 
the ISEE- A and ISEE-C spacecraft; the European Space Agency is responsible 
for ISEE-B. 

Each of the ISEE spacecraft, considered alone, can contribute to scientific 
knowledge. However the scientific basis for the mission relies on the spatial 
and temporal resolution obtained by comparing measurements made with iden- 
tical instruments on ISEE-A and ISEE-B. After launch, ISEE-C will supply 
information on the upstream solar wind, thereby enltancing the value of the 
total experiment. For maximum usefulness, this information should be ob- 
tained concurrently with measurements on ISEE-A and ISEE-B, A minimum 
of 2 years of simultaneous data is planned for the mission. 

ISEE-C Mission Orbit; An ecliptic plane projection of the nominal ISEE-C trans- 
fer trajectory is presented in Figure 1. Following a retro maneuver at approxi- 
mately injection plus 107 days, ISEE-C will be placed in a roughly elliptical 
path, termed a halo orbit, about the unstable Sun- Earth interior libration point 
. The nominal period of the halo orbit is approximately 6 months. 

Figure 2 depicts the halo orbit, and defines a Rotating Libration Point (RLP) 
coordinate system centered at . The X axis points toward the Earth-Moon 
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Figure 2. Halo Orbit ArouiKi #5un-Farth Libration Point 



barycenter, the Z axis points toward the north ecliptic pole (NEP), and the 
Y axis completes the right-handed coordinate frame. In the RLP, the maximum 
excursions in the position are X = +240,000 km, -170,000 km; Y = ±670,000 km; 
and Z = ±120,000 km. The stationkeeping strategy will be to maintain the space- 
craft within a torus of radius 5, 000 km about the nominal path. Stationkeeping 
maneuvers will probably be required every 43 to 60 days. 

ISEE-C Spacecraft; ISEE-C is spin stabilized, with a nominal spin rate of 
19.75 ±0.2 rpm. As shown in Figure 3, the spacecraft is a 16-slded, drum- 
shaped structure. It is 161 centimeters high, 174 centimeters in diameter, 
and has a nominal mass of 457 kilograms. Solar arrays cover the sides of 
the spacecraft except for a band near the center from which two experiment 
booms, t^vo Inertia booms, and four deployable wire antennas extend. These 
four wires combined with deployable ±Z axis wire antennas provide a three- 
dimensional radio mapping capability. 

The tolerance on the spacecraft dynamic balance is small. Prior to deploy- 
ment of the four radial wire^, the maximum spacecraft dynamic Imbalance 
is expected to be less than 0.2 degree. After deployment of the wires, the 
dynamic Imbalance is expected to decrease to approximattl: 0.1 degree. 

ISEE-C Attitude Sensors; A panoramic attitude scanner (PAS) and two fully 
redundant fine Sun sensors (FSSs) are located on the sides of the spacecraft. 

The PA3, manufactured by Ball Brothers Research Corporation, is Identical 
to the one flown on the International Ultraviolet Explorer and the ISEE-A 
spacecraft. The function of the PAS is to map the horizons and terminators 
of the Earth and Moon by measuring the crossing times relative to a Sun 
crossing. 

The FSS system is manufactured by the Adcole Corporation. It is similar to 
other Adcole 'Ugital Sun sensors, employing a quartz block with optical mask- 
ing of an array of silicon photocells and having a command slit to generate 
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Figure 3. ISEE-C Configuration 
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the command pulse. The fan-shaped field of view of each FSS is specified to 
be 128 ± 2 degrees In width. The FSSs meaf *re the angle bettveen the sunline 
and a reference axis fixed in the sensor. This measurement is readily converted 
to the Sun angle* 8 * bet^veen the spacecraft spin axis and the sunline. 

The FSSs operate on a 5 Volt system. The output is biased, quantized in 20 
millivolt counts, and telemetered to the ground. In the range -1 , 

the sy^em provides for a resolution of 0.003 degree. The measurement noise 
is specified to be less than O.OOS degree. When the spacecraft is spinning be- 
tween 5 and 25 revolutions per minute, the accuracy of the Sun angle measure- 
ment is as follows. 

Range of 8 Accuracy 

(Degrees) (Degrees) 


23-50; 

130-154 

iO.25 

50-87; 

93-130 

iO.l 

87-93 


xO. 03 


The increased accuracy in the S7 ^ ^ 93 degree range is due to calibration 

of the sensor which was performed by the manufacturer. A calibration table 
was generated for both FSSs at 0. 1-degree intervals over the specified Sun 
angle range. Each point in the table gives a correction term (63) to be added 
to the measured Sun ai^le. The root mean square (RMS) residual errors 
following application of the calibration data have been computed to be 0. 006 
and 0.004 degree for FSS 1 and 2, respectively. The ma.xlmum residual errors 
following application of the calibration data for FSS 1 and 2 have been computed 
to be 0.03 and 0.01 degree, respectively. A plot of the residual errors for 
FSS 2 is presented in Figure 4. 

The FSS mounting tolerance is very small. Tlte serisors are to be aligned with 
respect to the spacecraft body axes to within 0.005 degree. As a result, it is 
anticipated that the only significant Sun angle bias source will be the afore- 
mentioned spacecraft djmamlc imbalance of 0.2 degree or less. 
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Figure 4. FSS Measure. Errors Following Caliluration 


ISEE-C Attitude Determination and Control RequlremetiJu fc The ISEF-C attitude 
control requirement in its mission orbit is to maintain the spin axis to within 
1 degree of the North Ecliptic Pole (NEP). The attitude is to be determined to 
an accuracy of 1 degree. 

The spin axis attitude will be perturbed as a result of both orbit maneuvers and 
solar radiation pressure torque. Studies have shown that solar radiation torque 
will cause the spin axis to process along a small circle on the celestial sphere 
with a period of 1 year. The angular radius of the small circle at the start of 
the halo orbit is predicted to be approximately 0.3 c »gree. This results in a 
precession rate of approximately 0.008 degree per daji^. The angular radius of 
the precession circle is expected to decrease throughout the mission lifetime 
because of the shift in the center of mass relative to the center of pres ^e 
resulting from fuel usage. A precession radius of approximately 0.2 degree 
is anticipated when the tanks are empty. 

Because of the halo orbit geometry, the quality of the attitude information pro- 
vided by the PAS will vary significantly with position in orbit. In order to meet 
the 1 degree attitude determination accuracy requirement, up to 30 days of 
FSS data at an attitude perturbed only by the solar radiation pressure torque 
will be processed. The techniques described below that have been developed 
to perform this task are termed FSS-only techniques. 

Coordinate Systems; To perform the FSS-only computations, three coordinate 
frames are required: (1) the Geocentric Equatorial Inertial IBci) coordinute 
system, (2) the Geocentric Solar Ecliptic (GSE) coordinate system, and ( 2 ) the 
Geocentrlo Sun Llne-of-Slght Rotating (GSR^ coordinate system. 

The GCI coordinate system is Earth centers'’ The X_ axis points toward 

GCI 

the vernal equinox of date; the Z _ axis is normal to the true equator date; 

GCI 

and the Y... axis completes the right-handed system. The sunline and 



spacecraft vectors in GCI at time t are represented by U and H , respec> 
tively. (The unit sunline vector in GCI is U . ) The corresponding velocity 

vectors are U and 6 » respectively. Relative to the spacecraft, the sunline 
position and velocity vectors are 


u* ^ = u - 

A 

(1) 

A 

(2) 

The corresponding unit vectors are 


^A “ 1 U - R 1 

'i' U - R 

u* = -i 

(3) 

(4) 


The GSE w’oorJinate system is defined as follows. The X___ axis points 

Go£ 

toward the Sun. The axis points toward the north ecliptic pole (NEP). 

The Yggj, axis completes the right-handed system. The GCI coordinate 
system is related to the GSE coordinate system by the following transforma- 
tion matrix. 



sin L cos € 
cos L cos e 
-sin € 


sin L sin c 
cos L sin c 
cos e 



(5) 


where L is the celestial longitude of the Sun, and € is the obliquity of " e 
ecliptic. 



The GSR coordinate system is defined as follows. The axis is aligned 

/\ GSR 

with U (the geocentric unit vector parallel to the llne-of-sight from the 

A 

spacecraft to the Sun), The axis is defined by the following expression 


-A 

Y 


GSR 


A A 

NEP X U, 
A 

InSp X i> ! 


(6) 


where NEP is a unit vector pointing toward the north ecliptic pole. The 

Z_„ axis completes the right-handed system. 

GSR 

Referring to Figure a, the unit spin axis vector, A , is defined by the coor- 


dinates 8 and EL . 8 is the aspect angle of the sunline relative to A , EL 

is the dihedral angle between the X_„ - Z plane and the plane defined by 

GSR GSR 

A <A 

A . EL is positive in the sense Y___ x Z , 

GoR GSR GSR 


It follows that 


A = cos S Xjjjp - sin S sin EL + sin 3 cos EL 


( 7 ) 


The GSR coordinate system is aligned with the GSE coordinate system by using 
tlie following transformation matrix 



N 'E "E 

cos 


I sin y^ cos y^ 


- sin y 


N 


0 



( 8 ) 
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GSR Coordinate System 



where 


y = -sin"^ (U , • NEP) 
^ A 


(9) 




-1 


NEP X U 


I NEP XU 1 

A 


A ^ 

(NEP X U) 


( 10 ) 


FSS-Only Attitude Determination Equations; We denote time-dependence for 
an arbitrary vector, V , as V(t) . Figure 6 depicts the GSR coordinate frame 
on the unit celestial sphere at arbitrary time t . Here, the size of the locus 
of the sunline relative to the spacecraft is exaggerated for clarity. In reality, 
the maximum angular deviation of U . (t) above or below the ecliptic for the 

At 

ISEE-C mission orbit is approximately 0.05 degree. The maximum angular 

/\ .-s 

deviation of from the true sunline U(t) as measured in the ecliptic is 

approximately 0.25 degree. 

The normal to the instantaneous plane of motion of U . (t) is computed as 

A 

follows. 


n(t) = 


l-’A(t) X U^(t) 

I u^(t) X U^(t) i 


( 11 ) 


The angle C between the instantaneous plane of motion of U . (t) and the 

A 

^GSR “ 'gSR 


C - -sin Vn(t) . 


NEP X U (t) 


NEP X U , (t) 
A 


( 12 ) 
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The instantaneous angular velocity of l\(t) is 

A 


= cos"^|(U^(t) +'u^{t)) * V 


'a<') 


(13) 


Holding the GSR frame fixed at t , the unit suioline relative to the spacecraft 
at time t^ (where t^ = t -f- 6t) is computed as follows. 


l\(t) = cos + cos c sin 

+ Sin Z sin (w^^(t^-t)) 


(14) 


In the same frame of reference, the spin axis at time t is 


'^..l 


A(t ) = cos (8 (t) + 6^) X^gj^ - sin (8 (t) + 6^) sin (EL (t) + 6^) 

+ sin (8 (t) + 6^) cos (EL (t) + 6^) Z^g^^ (15) 


where 6^^ and 6^ are perturbations in S(t) and EL(t) , respectively, re- 
sulting from solar radiation pressure torque. 

E'sing Equations (14) and (15), the Sim angle at time t^ is computed to be 


5(tS = cos ^ [cos (d(t) + 6^) cos (o>^(t^-tn 

- sin (/3(t^ (EL(t) + 5^) cos f sin iuj^(t^-t)) 

+ sin (/?(t) - cos (EL(t^ - 6^) sin t; sin (a}^(t^-tn] 
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Differentiating with respect to t 


1 


cos ijSit) + 6 ) sin 

9 1 A i A 

sin (^(t) 6 ^) sin (EL(t) + 6.,) cos J cos ( u»^(t t)) (16i 

- sin ip(t) + 6j^) cos (EL(t> + 6^) sin £ cos U'] / sin/3(t^) 


Note tliat the perturbations fl and 6 are small relative to the motion of 

X M 

the sunline. By taking the limit of Equation (16) as 6t approaches zero, we 
obtain 


\ • 

lim — ,S(t) = o) cos Z sin EL(t) - to sin Ceos EL(t) (17) 
6t-0 3t^ ^ 


which becomes, after rearranging terms, 


,a(t) 

sin EL(t) = 7 + tan C cos EL(t) (18) 

’ 0 ) cos C 


Equation (IS) in general is solved as a quadratic in terms of sin EL^t). How- 
ever, as noted pre\iousIy, the ISEB-C spin axis will be aligned within 1 degree 
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of the NEP . Thus EL(t) is a small angle, and it follows that 


EL(t) 


U). COS 
A 


• + tan 


C 


(19) 


From the equation, it is evident that; (1) the instantaneous value of EL is 
directly proportional to the Instantaneous time rate of change of $ , and (2) 

Sun angle bias errors do not affect the computation of EL. 

Simulations have shown that the ISEE-C halo orbit Sun angle measurements 
over a full year can be accurately fit to a sine wave model Including the funda- 
mental and second harmonic frequencies. The linearized form of the analj'tical 
curve is as follows. 

/3(t) = Cq + cos wt + Cg sin wt + cos 2 wt + sin 2 tot (20) 

where to = 1 revolution per year and t is measured relative to an attitude 
reference time. Differentiating tvith respect to time yields 

$ (t) = -toCj sin tot + twCg cos tot - 2 toC^ sin 2 tot + 2 uJC^ cos 2 tot (21) 

Shorter arcs of data can be fit using either a sine wave model with fundamental 
frequency only, or using a low-order polynomial. Regardless of the analytical 
model employed, a smoothed valu : of i8(t) and 0{t) at arbitrary time t 
Nvithin the Interval of data can be evaluated. Attitude determination then pro- 
ceeds as follows. 

1. Compute EL(t) using Equations (12), (13), and (19) 

2, Compute A(t) in GSR using Equation (7) 
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3. 


to GSE using Equation (8) 


4. 


0 m 


Transforna A(t) 

\ 

Transform A(t)^ to GCI using the inverse transformation de- 
GoE 

fined by Equation (5) 

A(t) can then be converted to right ascension and declination 
GCI 

using the standard equations 


FSS-Only Attitude Uncertainty; Referring to Equation (19), the instantaneous 
elevation angle is a fimction of C , to , and 8 . Both C and to can be 

A A 

computed from the spacecraft ephemeris using Equations (12) and (13). Analy- 
sis has shown that the ISEE-C elevation angle computed in mission orbit will not 

be sensitive to eri-ors in either Z or to . . The maximum error in EL cor- 

A 

responding to the aitticipated uncertainties in spacecraft ephemeris will be less 
th:m 0. 01 degree. The uncertainty in EL will thus be a function of <3 uncer- 
tainty alone. The variance in EL is computed as follows. 

9 


8 


EL 2 2 ^ 

to. cos Z 
A 


( 22 ) 


2 - T- r 

where crj = Lh. ] [ P ] [h . J 
P 1 1 

[P] is the covariance matrix of error in the analytic curve coefficients 


[h.l - . where C. are the analytic cune coefficients 

11 i 

Note tliat [P] should be computed using the maximum relative uncertainty in 

8 measurements, cr^ This will avoid an unrealistically 

oME.Abl REAIENT 

large weighting of the Sun angle measurements. 

The variance in 3 can be computed as follows 


•> 





BIAS 


(23) 


2al 



where 


LgjJ = 3 ^/cC. 

0 

^ = Sun angle bias variance 

6 LAS 

The altitude uncertainty (standard deviation of the arc length between the true 
and estimated spin axis attitude vectors) is 

o’ 

<x“ -f a: (24) 

ARC EL ^ 

FSS-Only Prelaunch Attitude Uncertainty Predictions; ISEE-C FSS-only pre- 
launch attitude uncertainty predictions are presented in Figure 7 as a function 
of data span and data rate. Each Sun angle sample consists of a group of aver- 
aged measurements such that the effects of measurement noise and residual 
spacecraft nutation (following maneuvers) are minimized. A value of 
0. 03 degree was used for the relative uncertaintj’ in ^ measurements, 

Qo , in the computation of ihe covariance inatrLx of error lPj . 

^lEASUREMENT 

This corresponds to the maximum error obseiwed following sensor calibration. 
A value of 0.20 degree was used for • This corresponds to the afore- 

mentioned maximum anticipated offset betvveen the spacecraft principal axes 
and geometric axes. 

Referring to the figure, it is seen that the ISEE-C attitude can be determined 
to within an accuracy of 1 degree or better using data arcs of at least 10 days. 
.Attitude accuracy is a function of data rate for shorter data arcs. Improved 
accuracies are obtained usit^; the higher data rates. However, for data arcs 
exceeding 20 days, attitude accuracy cannot be improved using higher data 
rates. 

.Applicability of FSS-Only Technique to Other Missions; Equation (19) can be 
solved for 0 as follows. 

0 = oj. cos C sin EL - tan C cos EL (25) 

A 
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Differentiating with respect to EL , one obtains the following. 


cos C cos EL tan C sin EL 
Taking the Inverse of Uie above expression. 


1 

bp (u . cos C cos EL * tan C sin EL 
A 


t2G) 


(27) 


For ISEE-C, the motion of Uie spacecraft relative to the sunline is small. As 
a result, the maximum value of C Is on the order of 0. 1 degree. By inspec- 
tion of Equation i27), it is seen that small errors in will be transformed 
into large errors in the computation of EL when EL approaches 90 degrees. 
For example, when El. equals 60~degrees, the error magnification factor 
is approximately 2. Thus, the FSS-only attitude determination technique is 
not applicable to ISEE-C when EL is very large. 

There are other missions to which the FSS-only technique can be applied. For 
example, spin stabilized geosjmehronous satellites are often oriented toward 
the north celestial pole. -Mso, the solar radiation pressure torque effect is 
normally small. For such missions, EL will vary between ±23.5 degrees. 

The ma.ximum value of C will be less than 0.4 degree. The maximum 3 
error magnification factor computed using Equation (27) is approximately 1.1. 

It is concluded that the FSS-only attitude determination teclmlque can be applied 
to these missions. The attitude letermination accuracy obtained will, of 
course, be a function of the relative accuracy of the Sun sensor. 

Summary and Conclusions: Techniques have been developed to determine the 
ISEE-C spin axis attitude using a finite arc of smoothed Sen angle data, the 
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tedmiques differ from more conventional attitude determiiiation methods. in 
that ephemeris vectors are acquired onlj’ to transform the attitude results 
standard coordinate frames. The slowly- varying ISEE-C spin axis attitude 
dynamics resulting from solar radiation pressure toi'que are automatically 
accounted for in the attitude determination process. 

It is assumed that the ISEE-C FSS system will perform in orbit according to 
specifications. It is also assumed that the FSS mounting tolerances and 
ISEE-C spacecraft djaiamic Imbalance tolerances will not be exceeded. Given 
these assumptions, the ISEE-C spin axis attitude can be determined to within 
1.0 degree with the FSS-only techniques using up to 30 days of FSS data. The 
FSS-only techniques can also be applied to other missions, such as spin stabi- 
lized geosynchronous missions. 
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W79-14135 

INFRARED HORIZON SCANNER ATTITUDE 
DATA ERROR ANALYSIS FOR SEASAT-A* 

M. C. Phenneger, C, Manders, 

C. B. Spence, Jr., M. Levitas, 
and 

G. M, Lerner 
ABSTRACT 

This talk presents the results of a study of the effect of variations in the Earth’s 
seasonal and geographical horizon ladiance on the location of the infrared horizon 
as measured by ITHACO scanwhocls. Two types of variations are considered. 
These are (1) systematic variations of the mean (averaged over all longitudes) 
atniosphcric radiance duo to macroscopic changes in temperature as a function 
of latitude and season and (2) random variations in atomsphcric radiance due 
to microscopic fluctuations (weather). The effect of 'variations in the scanner 
wheel speed.- on the attitude determination accuracy is also presented. The 
computed horizon radiance and wlieel speed variation - induced attitude errors 
are than combined with errors caused by sensor alignment and electronics 
tolerances to obtain an overall estimate of the Scasat-A pitch and roll angle 
accuracy. 


♦ York supported by the Attitude Determination and Control Section, Goddard 
Space Flight Center , National Aeronautics and Space Admbiistration, imder 
Contract No. NAS 5-11999. 

** Scanwheel is registered trademark of ITHACO, Inc. 
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PATH OF HORIZON SCANNER ACROSS THE EARTH 
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NOTE SHADING INDICATES ROTATING FLYWHEEL AND PRISM ASSEMBLY 
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SEASAT OPTICS SPECTRAL RESPONSE 
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THE EARTH IR EMISSION FROM THE SAHARA AND THE ANTARCTIC 
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SIMULATION OF THE COLD CLOUD EFFECT 
ON THE EARTH'S RADIATION SPECTRUM 
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AVERAGE AND COLD CLOUD RADIANCE PROFILES AT THE EQUATOR IN JULY 



TANGENT HEIGHT (KM> 









JULY PITCH ERRORS WITH COLD CLOU 
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ATTITUDE ACQUISITION CONTINGENCY STUDIES 

FOR 

THE APPLICATIONS EXPLORER MISSIONS- A/HEAT CAPACITY 
MAPPING MISSION (AEM-A/RCMM) SPACECRAFT 

\Vhittak Huang, Mihaly G. Grell 
and 

Gerald M. Lemer 
Compxiter Sciences Corporation 

ABSTRACT 

H)C Heat Capacity Mapping Mission (ilCMM) is the first of a series of satellites 

utilizing a basic, modularly designed launch vehicle and satellite support system 

called the Applicatiems Exiilorer Missions (AEM) HCMH to be 

a 

launched in April 1978 into a 600 km altitude. Sun- synchronous polor orbit, will 
conduct a thermal mapping of the Noith American continent to investigate Earth 
resources' availability. The spacecraft has an attitude control system consisting 
of wheel-mounted infrared horizon sensor oriented along the negative body Y-axis 
<thc orbit normal for the nominal attitude), a 3-axis magnetometer and 3-orUiogonal 
electromagnetic coils. The magnetometer data is used for mission-mode 3-axis 
attitude control (pitch = roll= yaw = 0). Control laws, first proposed by Saymor 
Kant, Peter llui and Joseph Lidston , which relate the attitude data from 

the sensors to the control torque commands arc used by the onboard attitude 
computer to achieve attitude acquisition and to maintain the mission-mode 
attitude. 

Attitude acquisition requires maneuvering the spacecraft from the attitude after 
separation (spinning about the body Z-axis aloi^ the velocity vector) to the mission 
mode attitude. Attitude maintenance consists of pitch control by modulating the 
scanner wheel speed, momentum control by commanding the X- and Z-axis 
elcctromagcnts, and roll and r.utation control by commanding the V-axis electro- 
magnet. Yaw is controlled indirectly via gyrocompassing. 
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